Tweet

これは Math Advent Calendar 2017 6日目の記事です。

はじめに

ベイズBPM計測器は、曲に合わせて手動でスペースキーを押すとBPMを測定してくれるツールです。

今日は、スペースキーが押された時刻から、ベイズの定理を用いてBPMを推定する方法について考えたいと思います。

観測データ

まず、スペースキーが入力された時刻を

$$T_n \quad (0 \leq n \leq N)$$

（ミリ秒）とし、スペースキーが入力された時刻の間隔を

$$a_n = T_{n+1} - T_n \quad (0 \leq n \leq N - 1)$$

と定義します。

簡単な推定方法

BPMの簡単な推定方法として、最後にキーが押された時刻と最初にキーが押された時刻の差を取って、(キーが押された回数 - 1) で割る、という方法があります。

$$E = \frac{1}{N} \times ( T_N - T_0 ) = \frac{1}{N} \times \sum_{i=0}^{N-1} a_i$$

しかしこの方法では、最初のキー入力と最後のキー入力以外の時刻が計算結果に影響していません。これは得られた情報を最大限に利用してBPM推定するアルゴリズムとは考えづらいと予想できます。本当でしょうか？ では考えてみましょう。

事前分布

楽曲の本来のBPMから計算される1拍の長さを $2u$ 、最初の拍と最後の拍が鳴る時刻の（相加）平均を $t$ とすると、$i$ 拍目の時刻は

$$t + (-N + 2i)u \quad (0 \leq i \leq N)$$

となります。ここで、t と u が事前分布として一様分布を取る場合を考えます。つまり、BPM 180の曲が与えられる確率と、BPM 150の曲が与えられる確率と、BPM 120の曲が与えられる確率がどれも等しいということを意味します。これは、ある程度妥当な仮定だと考えられます。

スペースキーを人間が押す時刻が、この時刻を中心とする独立で分散が $\sigma^2$ の正規分布に従うと仮定すると、時刻 $x$ に $i$ 回目のキーが入力される確率は次のように表されます。

$$p_i(x) = \mathcal{N}(t - Nu + 2iu, \sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left( -\frac{\{x - (t - Nu + 2iu)\}^2}{2 \sigma^2} \right)$$

N=3 の場合

話を簡単にするために、N = 3 と固定して考えてみましょう。この場合、0、1、2、3拍目の時刻はそれぞれ

$$t - 3u, \quad t - u, \quad t + u, \quad t + 3u$$

となります。従って、時刻 $x$ に0、1、2、3回目のキーが入力される確率はそれぞれ次のように表されます。

$$\begin{split} p_0(x) &= \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left( -\frac{\{x - (t - 3u)\}^2}{2 \sigma^2} \right) \\ p_1(x) &= \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left( -\frac{\{x - (t - u)\}^2}{2 \sigma^2} \right) \\ p_2(x) &= \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left( -\frac{\{x - (t + u)\}^2}{2 \sigma^2} \right) \\ p_3(x) &= \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left( -\frac{\{x - (t + 3u)\}^2}{2 \sigma^2} \right) \end{split}$$

また、実際に観測されたキー入力の時刻は、次のように表されます。

$$T_0, \quad T_1, \quad T_2, \quad T_3$$

ここで、一般性を失わずに $T_0 = 0$ とすると、スペースキーが入力された時刻の間隔 $a_i$ を用いて次のように表すことができます。

$$\begin{split} T_0 &= 0 \\ T_1 &= a_0 \\ T_2 &= a_0 + a_1 \\ T_3 &= a_0 + a_1 + a_2 \end{split}$$

ある t と u の値が与えられた場合にこの実現値が観測される確率密度関数は

$$p(a_0, a_1, a_2 \, | \, t, u) = p_0(0) \cdot p_1(a_0) \cdot p_2(a_0 + a_1) \cdot p_3(a_0 + a_1 + a_2)$$

と表すことができます。

事後確率と1拍の長さの計算

t と u を、閉区間 [-L, L] で一様な値を取る分布とすると、このことは次のような式で表すことが出来ます。

$$p(t, u) = \frac{1}{4L^2} \quad (-L \leq t, u \leq L)$$

また、この範囲の任意の t, u に対して $a_0, a_1, a_2$ が観測される確率は、次のように計算できます。

$$p(a_0, a_1, a_2) = \int_{-L}^{L} \int_{-L}^{L} p(a_0, a_1, a_2 \, | \, t, u) \, p(t, u) \, \mathrm{d}t\mathrm{d}u$$

さて、ベイズの定理より、

$$p(t, u \, | \, a_0, a_1, a_2) = \frac{p(a_0, a_1, a_2 \, | \, t, u) \, p(t, u)}{p(a_0, a_1, a_2)}$$

となりますので、上の式を代入して約分すると次のようになります

$$p(t, u \, | \, a_0, a_1, a_2) = \frac{p_0(0) \cdot p_1(a_0) \cdot p_2(a_0 + a_1) \cdot p_3(a_0 + a_1 + a_2)}{\int_{-L}^{L} \int_{-L}^{L} p_0(0) \cdot p_1(a_0) \cdot p_2(a_0 + a_1) \cdot p_3(a_0 + a_1 + a_2) \, \mathrm{d}t\mathrm{d}u}$$

ところで、この確率密度関数は $a_0, a_1, a_2$ が観測された場合の t, u の確率密度関数ですから、u に関して期待値を取れば、1拍の長さの推定値が得られることになります。u の期待値を計算すると、次のようになります。

$$E = \int_{-L}^{L} \int_{-L}^{L} p(t, u \, | \, a_0, a_1, a_2) \cdot u \, \mathrm{d}t\mathrm{d}u$$

実際に計算してみた

では、まずは分母の $p(a_0, a_1, a_2)$ から計算していきましょう。被積分関数は $p(a_0, a_1, a_2 \, | \, t, u)$ と $p(t, u)$ の積ですが、前者は

$$\begin{split} p(a_0, a_1, a_2 \, | \, t, u) &= p_0(0) \cdot p_1(a_0) \cdot p_2(a_0 + a_1) \cdot p_3(a_0 + a_1 + a_2) \\ &= \frac{1}{(2\pi\sigma^2)^2} \exp\left( -\frac{\{0 - (t - 3u)\}^2}{2 \sigma^2} \right) \exp\left( -\frac{\{a_0 - (t - u)\}^2}{2 \sigma^2} \right) \\ & \quad\quad\quad\quad \times \exp\left( -\frac{\{a_0 + a_1 - (t + u)\}^2}{2 \sigma^2} \right) \exp\left( -\frac{\{a_0 + a_1 + a_2 - (t + 3u)\}^2}{2 \sigma^2} \right) \\ &= \frac{1}{(2\pi\sigma^2)^2} \exp\left( -\frac{\{0 - (t - 3u)\}^2 + \{a_0 - (t - u)\}^2 + \{a_0 + a_1 - (t + u)\}^2 + \{a_0 + a_1 + a_2 - (t + 3u)\}^2}{2 \sigma^2} \right) \\ &= \frac{1}{(2\pi\sigma^2)^2} \exp\left( -\frac{ 3a_0^2 + 2a_1^2 + a_2^2+ 2(a_0a_1 + a_1a_2 + a_2a_0) + 4t^2 - 2(3a_0 + 2a_1 + a_2)t + 20u^2 - 2(3a_0 + 4a_1 + 3a_2)u }{2 \sigma^2} \right) \end{split}$$

となるので、ゴリゴリと計算すると次のようになります。（一部で途中式を省略しています）

$$\begin{split} & p(a_0, a_1, a_2) \\ &= \int_{-L}^{L} \int_{-L}^{L} p(a_0, a_1, a_2 \, | \, t, u) \, p(t, u) \, \mathrm{d}t\mathrm{d}u \\ &= \frac{1}{4L^2} \int_{-L}^{L} \int_{-L}^{L} p(a_0, a_1, a_2 \, | \, t, u) \, \mathrm{d}t\mathrm{d}u \\ &= \frac{1}{4L^2 (2\pi\sigma^2)^2} \int_{-L}^{L} \biggl\{ \exp\left(-\frac{ - \frac{(3a_0 + 2a_1 + a_2)^2}{4} + 3a_0^2 + 2a_1^2 + a_2^2+ 2(a_0a_1 + a_1a_2 + a_2a_0) + 20u^2 - 2(3a_0 + 4a_1 + 3a_2)u }{2 \sigma^2}\right) \\ & \quad\quad \times \int_{-L}^{L} \exp\left( -\frac{4}{2 \sigma^2}\left(t - \frac{3a_0 + 2a_1 + a_2}{4}\right)^2 \right) \, \mathrm{d}t \biggr\} \mathrm{d}u \\ &= \frac{1}{4L^2 (2\pi\sigma^2)^2} \\ & \quad\quad \times \int_{-L}^{L} \exp\left(-\frac{ - (3a_0 + 2a_1 + a_2)^2 + 4\{\cdots + 20u^2 - 2(3a_0 + 4a_1 + 3a_2)u\} }{8 \sigma^2}\right)\mathrm{d}u \\ & \quad\quad \times \int_{-L}^{L} \exp\left( -\frac{4}{2 \sigma^2}\left(t - \frac{3a_0 + 2a_1 + a_2}{4}\right)^2 \right) \, \mathrm{d}t \\ &= \frac{1}{4L^2 (2\pi\sigma^2)^2} \\ & \quad\quad \times \exp\left(-\frac{ - (3a_0 + 4a_1 + 3a_2)^2 - 5(3a_0 + 2a_1 + a_2)^2 + 20(3a_0^2 + 2a_1^2 + a_2^2) + 40(a_0a_1 + a_1a_2 + a_2a_0) }{40 \sigma^2}\right) \\ & \quad\quad \times \int_{-L}^{L} \exp\left( -\frac{80}{8 \sigma^2} \left(u - \frac{3a_0+4a_1+3a_2}{20} \right)^2 \right) \, \mathrm{d}u \\ & \quad\quad \times \int_{-L}^{L} \exp\left( -\frac{4}{2 \sigma^2}\left(t - \frac{3a_0 + 2a_1 + a_2}{4}\right)^2 \right) \, \mathrm{d}t \\ &= \frac{C}{4L^2 (2\pi\sigma^2)^2} \\\end{split}$$

ただしCは定数となります。ここではまだCは計算しません。

期待値の計算

次に $E$ を計算します。上の計算とほぼ一緒なので簡単です。（一部で途中式を省略しています）

$$\begin{split} E &= \int_{-L}^{L} \int_{-L}^{L} p(t, u \, | \, a_0, a_1, a_2) \cdot u \, \mathrm{d}t\mathrm{d}u \\ &= \int_{-L}^{L} \int_{-L}^{L} \frac{p(a_0, a_1, a_2 \, | \, t, u) \, p(t, u)}{p(a_0, a_1, a_2)} \cdot u \, \mathrm{d}t\mathrm{d}u \\ &= \frac{1}{4L^2} \cdot \frac{4L^2 (2\pi\sigma^2)^2}{C} \int_{-L}^{L} \int_{-L}^{L} p(a_0, a_1, a_2 \, | \, t, u) \cdot u \, \mathrm{d}t\mathrm{d}u \\ &= \frac{(2\pi\sigma^2)^2}{C} \int_{-L}^{L} \int_{-L}^{L} p(a_0, a_1, a_2 \, | \, t, u) \cdot u \, \mathrm{d}t\mathrm{d}u \\ &= \frac{1}{C} \int_{-L}^{L} \int_{-L}^{L} \exp\left( -\frac{ \cdots + 4t^2 - 2(3a_0 + 2a_1 + a_2)t + 20u^2 - 2(3a_0 + 4a_1 + 3a_2)u }{2 \sigma^2} \right) \cdot u \, \mathrm{d}t\mathrm{d}u \\ &= \frac{1}{C} \int_{-L}^{L} \biggl\{ \exp\left(-\frac{ - \frac{(3a_0 + 2a_1 + a_2)^2}{4} + \cdots + 20u^2 - 2(3a_0 + 4a_1 + 3a_2)u }{2 \sigma^2}\right) \\ & \quad\quad \times \int_{-L}^{L} \exp\left( -\frac{4}{2 \sigma^2}\left(t - \frac{3a_0 + 2a_1 + a_2}{4}\right)^2 \right) \, \mathrm{d}t \biggr\} \cdot u \, \mathrm{d}u \\ &= \frac{1}{C} \int_{-L}^{L} \exp\left(-\frac{ - (3a_0 + 2a_1 + a_2)^2 + 4\{\cdots + 20u^2 - 2(3a_0 + 4a_1 + 3a_2)u\} }{8 \sigma^2}\right) \cdot u \, \mathrm{d}u \\ & \quad\quad \times \int_{-L}^{L} \exp\left( -\frac{4}{2 \sigma^2}\left(t - \frac{3a_0 + 2a_1 + a_2}{4}\right)^2 \right) \, \mathrm{d}t \\ &= \frac{1}{C} \times \exp\left(-\frac{ - (3a_0 + 4a_1 + 3a_2)^2 - 5(3a_0 + 2a_1 + a_2)^2 + 20(3a_0^2 + 2a_1^2 + a_2^2) + 40(a_0a_1 + a_1a_2 + a_2a_0) }{40 \sigma^2}\right) \\ & \quad\quad \times \int_{-L}^{L} \exp\left( -\frac{80}{8 \sigma^2} \left(u - \frac{3a_0+4a_1+3a_2}{20} \right)^2 \right) \cdot u \, \mathrm{d}u \\ & \quad\quad \times \int_{-L}^{L} \exp\left( -\frac{4}{2 \sigma^2}\left(t - \frac{3a_0 + 2a_1 + a_2}{4}\right)^2 \right) \, \mathrm{d}t \\ &= \frac{ \int_{-L}^{L} \exp\left( -\frac{10}{\sigma^2} \left(u - \frac{3a_0+4a_1+3a_2}{20} \right)^2 \right) \cdot u \, \mathrm{d}u }{ \int_{-L}^{L} \exp\left( -\frac{10}{\sigma^2} \left(u - \frac{3a_0+4a_1+3a_2}{20} \right)^2 \right) \mathrm{d}u } \\ &= \frac{ \int_{-L}^{L} \frac{\sqrt{20}}{\sqrt{2\pi\sigma^2}} \exp\left( -\frac{20}{2\sigma^2} \left(u - \frac{3a_0+4a_1+3a_2}{20} \right)^2 \right) \cdot u \, \mathrm{d}u }{ \frac{\sqrt{20}}{\sqrt{2\pi\sigma^2}} \int_{-L}^{L} \exp\left( -\frac{20}{2\sigma^2} \left(u - \frac{3a_0+4a_1+3a_2}{20} \right)^2 \right) \mathrm{d}u } \\\end{split}$$

このまま L を $\infty$ に近づけると次のようになります。

$$E = \frac{ \int_{-\infty}^{\infty} \frac{\sqrt{20}}{\sqrt{2\pi\sigma^2}} \exp\left( -\frac{20}{2\sigma^2} \left(u - \frac{3a_0+4a_1+3a_2}{20} \right)^2 \right) \cdot u \, \mathrm{d}u }{ \frac{\sqrt{20}}{\sqrt{2\pi\sigma^2}} \int_{-\infty}^{\infty} \exp\left( -\frac{20}{2\sigma^2} \left(u - \frac{3a_0+4a_1+3a_2}{20} \right)^2 \right) \mathrm{d}u } \\$$

ここで分子は正規分布の期待値なので、簡単に求められます。

$$\int_{-\infty}^{\infty} \frac{\sqrt{20}}{\sqrt{2\pi\sigma^2}} \exp\left( -\frac{20}{2\sigma^2} \left(u - \frac{3a_0+4a_1+3a_2}{20} \right)^2 \right) \cdot u \, \mathrm{d}u = \frac{3a_0+4a_1+3a_2}{20}$$

また、分母はガウス関数の積分公式より次のようになります。

$$\int_{-\infty}^{\infty} \exp\left( -\frac{20}{2\sigma^2} \left(u - \frac{3a_0+4a_1+3a_2}{20} \right)^2 \right) \mathrm{d}u = \frac{\sqrt{2\pi\sigma^2}}{\sqrt{20}}$$

従って、これらを代入すると

$$E = \frac{3a_0+4a_1+3a_2}{20}$$

となります。

計算結果の考察

1拍の長さを 2u と仮定しましたので、求める1拍の長さの期待値はこの2倍となります。

$$2E = \frac{3a_0+4a_1+3a_2}{10}$$

これは、1拍目と2拍目の間隔と、2拍目と3拍目の間隔と、3拍目と4拍目の間隔を、3:4:3で混合した値となっていることがわかります。

Nの値が一般の場合

Nの値が一般の場合は、同様に計算すると、次のようになります。

$$2E = \frac{6}{N(N+1)(N+2)} \times \sum_{i=0}^{N} (2i-N)T_i$$

$$\mathrm{BPM} = \frac{60000}{2E}$$

または、同値ですが、次のように表すこともできます。

$$2E = \frac{6}{N(N+1)(N+2)} \times \sum_{i=0}^{N-1} (N-i)(i+1)a_i$$

$a_i$ の係数の具体的な値は次のようになります。

完成品

できたものはこちらに置いてあります。

[数学] ベイズBPM計測器を公開しました | yuinore.​net https://t.co/81EDvpWgik

記事を投稿しました。

— ゆいのあ (@yuinore) 2017年10月16日

明日は y_e_af さんで、「関数環P(X), R(X), A(X)について」です。お楽しみに。

Tweet