これは Math Advent Calendar 2017 6日目の記事です。
はじめに
ベイズBPM計測器は、曲に合わせて手動でスペースキーを押すとBPMを測定してくれるツールです。
今日は、スペースキーが押された時刻から、ベイズの定理を用いてBPMを推定する方法について考えたいと思います。
観測データ
まず、スペースキーが入力された時刻を
(ミリ秒)とし、スペースキーが入力された時刻の間隔を
と定義します。
簡単な推定方法
BPMの簡単な推定方法として、最後にキーが押された時刻と最初にキーが押された時刻の差を取って、(キーが押された回数 - 1) で割る、という方法があります。
しかしこの方法では、最初のキー入力と最後のキー入力以外の時刻が計算結果に影響していません。これは得られた情報を最大限に利用してBPM推定するアルゴリズムとは考えづらいと予想できます。本当でしょうか? では考えてみましょう。
事前分布
楽曲の本来のBPMから計算される1拍の長さを 、最初の拍と最後の拍が鳴る時刻の(相加)平均を
とすると、
拍目の時刻は
となります。ここで、t と u が事前分布として一様分布を取る場合を考えます。つまり、BPM 180の曲が与えられる確率と、BPM 150の曲が与えられる確率と、BPM 120の曲が与えられる確率がどれも等しいということを意味します。これは、ある程度妥当な仮定だと考えられます。
スペースキーを人間が押す時刻が、この時刻を中心とする独立で分散が の正規分布に従うと仮定すると、時刻
に
回目のキーが入力される確率は次のように表されます。
N=3 の場合
話を簡単にするために、N = 3 と固定して考えてみましょう。この場合、0、1、2、3拍目の時刻はそれぞれ
となります。従って、時刻 に0、1、2、3回目のキーが入力される確率はそれぞれ次のように表されます。
また、実際に観測されたキー入力の時刻は、次のように表されます。
ここで、一般性を失わずに とすると、スペースキーが入力された時刻の間隔
を用いて次のように表すことができます。
ある t と u の値が与えられた場合にこの実現値が観測される確率密度関数は
と表すことができます。
事後確率と1拍の長さの計算
t と u を、閉区間 [-L, L] で一様な値を取る分布とすると、このことは次のような式で表すことが出来ます。
また、この範囲の任意の t, u に対して が観測される確率は、次のように計算できます。
さて、ベイズの定理より、
となりますので、上の式を代入して約分すると次のようになります
ところで、この確率密度関数は が観測された場合の t, u の確率密度関数ですから、u に関して期待値を取れば、1拍の長さの推定値が得られることになります。u の期待値を計算すると、次のようになります。
実際に計算してみた
では、まずは分母の から計算していきましょう。被積分関数は
と
の積ですが、前者は
となるので、ゴリゴリと計算すると次のようになります。(一部で途中式を省略しています)
ただしCは定数となります。ここではまだCは計算しません。
期待値の計算
次に を計算します。上の計算とほぼ一緒なので簡単です。(一部で途中式を省略しています)
このまま L を に近づけると次のようになります。
ここで分子は正規分布の期待値なので、簡単に求められます。
また、分母はガウス関数の積分公式より次のようになります。
従って、これらを代入すると
となります。
計算結果の考察
1拍の長さを 2u と仮定しましたので、求める1拍の長さの期待値はこの2倍となります。
これは、1拍目と2拍目の間隔と、2拍目と3拍目の間隔と、3拍目と4拍目の間隔を、3:4:3で混合した値となっていることがわかります。
Nの値が一般の場合
Nの値が一般の場合は、同様に計算すると、次のようになります。
または、同値ですが、次のように表すこともできます。
の係数の具体的な値は次のようになります。
i=0 | i=1 | i=2 | i=3 | i=4 | i=5 | |
---|---|---|---|---|---|---|
N=1 | 1 | |||||
N=2 | 2 | 2 | ||||
N=3 | 3 | 4 | 3 | |||
N=4 | 4 | 6 | 6 | 4 | ||
N=5 | 5 | 8 | 9 | 8 | 5 | |
N=6 | 6 | 10 | 12 | 12 | 10 | 6 |
完成品
できたものはこちらに置いてあります。
[数学] ベイズBPM計測器を公開しました | yuinore.net https://t.co/81EDvpWgik
記事を投稿しました。
— ゆいのあ (@yuinore) 2017年10月16日
明日は y_e_af さんで、「関数環P(X), R(X), A(X)について」です。お楽しみに。