ダーツにおける狙い目を検証してみた！【第二弾】(フーリエ変換を利用したアルゴリズム改良)

#はじめに
こんにちは。
以前のブログでダーツ( 01, CRICKET )の期待値計算を行ったのですが、フーリエ変換を利用することで計算の高速化&計算精度の向上が実現したので、手法等を紹介したいと思います。
計算精度の向上により、 01 において最適な狙い目が BULL から TRIPLE20 へと切り替わる転移点が新たに確認できました。

👇前回のブログはこちら👇

https://qiita.com/takumisk/items/b1966c33aa5155c0f2c7

#ダーツのルール,期待値計算の準備
今回のブログでは 01 のみを考えます。得点に絞って解説すれば下図の通りです。
![Darts_value_01.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/713550/aca4eed0-42fe-6b94-5864-8206e436529e.png)
(前回のブログより引用)

以下ではブルの中心を原点に取る座標系を考え、上図を2次元平面上の関数 $f(x,y)$ として定義します(ボードの外は 0 点です)。

次に、常に狙った位置に刺さるとは限らないので、狙い目から $(x',y')$ だけズレた位置に刺さる確率を $p(x',y')$ とします。前回のブログに引き続き、以下では確率分布は2次元正規分布に従うとします。

```math
確率密度関数:p(x',y')=\frac{1}{2\pi\sigma^2}\exp\left(-\frac{r'^2}{2\sigma^2}\right) \quad\left(r'=\sqrt{x’^2+y'^2}\right)
```

$r'$　は狙うポイントからの距離、パラメータ　$\sigma$　はブレの度合いをそれぞれ表します。

以上の設定の下で、ある点$(x,y)$を狙った時の点数の期待値 $E(x,y)$ は以下のように計算されます。

```math
E(x,y)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(x',y')p(x-x',y-y') dx'dy'
```

数学的には、このような積分は関数 $f(x,y),p(x,y)$ の畳み込みと呼ばれます。畳み込み演算はフーリエ変換と相性がよいことが知られているので、以下にこの方針を検討します。

#数学的準備
##重積分の離散化
積分を離散化する際の隣接点同士の間隔を $d$ とします。

```math
f[n_x][n_y]=f(n_xd,n_yd), p[n_x][n_y]=p(n_xd,n_yd)d^2
```

と定めれば、近似的に $E[n_x][n_y]=E(n_xd,n_yd)$ は以下のように計算可能です。

```math
E[n_x][n_y]\simeq\sum_{m_x=-\infty}^{\infty}\sum_{m_y=-\infty}^{\infty} f[m_x][m_y]p[n_x-m_x][n_y-m_y]
```

これは数列に関する畳み込みとなっています。

##巡回畳み込みへの式変形
さて、$f,p$ はいずれも遠方では　0　か非常に小さい値を取るので、近似的に有限項で打ち切ることが可能です。以下では、$f$ は $-N\leq n_x \leq N,-N\leq n_y \leq N$ 、$p$ は $-M\leq n_x \leq M,-M\leq n_y \leq M$ でのみ非ゼロの値を取りうるものとします。
さらに都合上、以下のように $g,q$ を定めます。

```math
g[n_x][n_y]=f[n_x-N][n_y-N],\,q[n_x][n_y]=p[n_x-M][n_y-M]
```

$g$ は $0\leq n_x \leq 2N,0\leq n_y \leq 2N$ 、$q$ は $0\leq n_x \leq 2M,0\leq n_y \leq 2M$ でのみ非ゼロの値を取りうることになります。

```math
D[n_x][n_y]=\sum_{m_x=-\infty}^{\infty}\sum_{m_y=-\infty}^{\infty} g[m_x][m_y]q[n_x-m_x][n_y-m_y]\,(g,q\,の畳み込み)
```

とすれば、${E[n_x][n_y]=D[n_x+N+M][n_y+N+M]}$ が成立します。また、$D$ は $0\leq n_x \leq 2(N+M), 0\leq n_y \leq 2(N+M)$ の範囲でのみ非ゼロの値をとりうることとなります。

唐突ですが、ここで $L\geq 2(N+M)$ なる $L$ を用いて、$0\leq n_x \leq L,0\leq n_y \leq L$ 上で $\tilde{g},\tilde{q}$ を周期 $L$ の関数として、次のように定義します。

```math
\tilde{g}[n_x][n_y]=
\begin{cases}
g[n_x][n_y]\,(0\leq n_x \leq 2N+1,0\leq n_y \leq 2N+1)\\
0\quad\quad\quad(それ以外)
\end{cases}
\\
\tilde{p}[n_x][n_y]=
\begin{cases}
q[n_x][n_y]\,(0\leq n_x \leq 2M+1,0\leq n_y \leq 2M+1)\\
0\quad\quad\quad(それ以外)
\end{cases}
```

このような操作はゼロパディング(0埋め)と呼ばれます。

```math
\tilde{D}[n_x][n_y]=\sum_{m_x=0}^L\sum_{m_y=0}^L \tilde{g}[m_x][m_y]\tilde{q}[n_x-m_x][n_y-m_y]\,(\tilde{g},\tilde{q}\,の巡回畳み込み)
```

とすれば、 $D,\tilde{D}$ は $0\leq n_x \leq 2(N+M), 0\leq n_y \leq 2(N+M)$ の範囲で一致するので、$\tilde{D}$ を考えれば良いわけです。 $\tilde{D}$ は有限和なので愚直に計算することも可能ですが、 全ての $(n_x,n_y)$ について計算する場合、計算量は $O(L^4)$　となります。

## 巡回畳み込みと離散フーリエ変換(DFT)
まず、離散フーリエ変換(DFT)の定義を確認します。簡単のため、一次元に限って考えます(多次元でも基本は同様)。一般の数列 $f[n]\,(0\leq n \leq N-1)$ に対して、DFT後の数列 $F[k]\,(0\leq k \leq N-1)$ を次のように定めます。

```math
F[k]=\sum_{n=0}^{N-1} f[n]\exp(-2\pi jkn/N)
```
また、離散フーリエ逆変換を用いることで $F[k]$ から $f[n]$ を復元すること
が可能です。

```math
f[n]=\frac{1}{N} \sum_{k=0}^{N-1} F[k]\exp(2\pi jkn/N)
```

さて、一般の数列 $f,g$ の巡回畳み込み $f\ast g$ を考えてみましょう。

```math
{\rm DFT}[f\ast g]=\sum_{n=0}^{N-1} \sum_{m=0}^{N-1}f[m]g[n-m]\exp(-2\pi jkn/N)\\
=\sum_{m=0}^{N-1} f[m] \exp(-2\pi jkm/N)\sum_{n=0}^{N-1} g[n-m]\exp(-2\pi jk(n-m)/N)\\
=F[k]G[k]\,(f,g は周期 N の関数として拡張している)
```

より、畳み込みのDFTが $f,g$ のDFTの積として表せることがわかりました。
ところで、長さ $N$ の数列のDFTやその逆変換は、高速フーリエ変換(FFT)というアルゴリズムにより、$O(N\log N)$で変換可能です。畳み込みを愚直に計算する場合、計算量は $O(N^2)$ ですが、$f,g$それぞれをDFTし、掛け算してから逆変換すれば$O(N\log N)$での計算が実現します。
今回求めるべき $\tilde{D}[n_x][n_y]$ は

```math
\tilde{D}[n_x][n_y]={\rm DFT}^{-1}[{\rm DFT}[\tilde{g}[n_x][n_y]]{\rm DFT}[\tilde{q}[n_x][n_y]]
```
と辿ることで、$O(L^2 \log L)$　での計算が可能となります。

#前回手法との比較
積分計算の離散化幅を $d$ とすると、前回手法では計算量は $O(d^{-4})$ でしたが、今回の手法では $O(d^{-2} \log(1/d))$ となっています。計算量の制約により、前回手法では $d = 5 \,{\rm mm}$　が限界でしたが、今回の手法では $d = 0.5 \,{\rm mm}$ にまで下げることができました。計算精度を保つためには $d$ を誤差のパラメータ $\sigma$ より数倍小さく保つ必要があり、前回手法では $\sigma< 15\,{\rm mm}$ の範囲を精度よく計算することが叶いませんでしたが、今回の手法ではこの範囲でも精度を担保して計算することができました。

#計算結果
まずは、01における最適戦略のσ依存性を確認します。
![Darts_01_table.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/713550/2338bc25-e4b2-ddd8-42c6-87273710e5b5.png)
前回の結果と相違はありませんが、新たに $\sigma<15\,{\rm mm} $ の計算結果が加わりました。この図から、$\sigma=10\,{\rm mm}$,　$\sigma=8\,{\rm mm}$　の間で最適な狙い目が BULL から TRIPLE20 に切り替わることがわかります。
stats(3投の合計点数)のσ依存性をグラフにすると下図の通りです。
![Darts_sigma_stats_theo_narrow.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/713550/9887fe5c-6991-981f-6021-bcd47bf9d7d6.png)
最適戦略が切り替わる点での　stats　は　140　を超えていますが、プロでも　stats　が　140　を超す方はほとんどいないようなので、　TRIPLE20　が最適戦略になることはほぼあり得ないと言えます。
しかしながら、プロの試合などではアウターブル(　BULL　のうち、外側の部分)を　25　点と数える方式もあります。この場合は様子が一変します。最適戦略のσ依存性は以下の通りです。
![Darts_01_table_pro.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/713550/c53d5004-6ea1-e441-d197-35dceac93efa.png)
この図から、$\sigma=19\,{\rm mm}$,　$\sigma=18\,{\rm mm}$　の間で最適な狙い目が TRIPLE19 から TRIPLE20 に切り替わることがわかります。いずれにせよ、　BULL　が最適戦略となることはないようです。$\sigma=19\,{\rm mm}$　はAフライト, stats 90 相当ですので、Aフライトでもアウターブル 25 点方式では TRIPLE19 が最適戦略になりうることとなります。これは意外な結果かなと思います。

#さいごに
今回のブログでは、フーリエ変換を利用した計算量改善の一例を取り上げてみました。
FFT(高速フーリエ変換)は　numpy　パッケージの関数を利用するだけでしたが、それでも実装はそれなりに面倒でした。実装コストに見合う有意義な計算結果が得られたかは微妙ですが、フーリエ変換によるアクロバティックな手法で愚直法の計算結果とぴったり合う様子は感慨深いものがありました。
最後までお読みいただきありがとうございました。