# "初期Seed厳選"における最適戦略の考察(2)
## 序文

前回は"初期Seed厳選"における最適戦略を考察するために, 目的Seedを得るまでにかかる時間の期待値を計算した.

ただし, あくまで期待値は期待値であって, 実際に目的Seedを得るまでの時間はバラつきが生じる. あまりにもバラつきが大きい戦略を取ってしまうと, いわゆる確率の"下振れ"を引いてしまったときに何時までたっても厳選が終わらない, といったことも考えられる.

今回は, このバラつきを考慮に入れて, 確率の下振れを引いた場合でもなるべく少ない時間で厳選を終わらせる方法を考察する.

## 事前条件
乱数生成器の周期が十分に大きく, 各出力が独立であると見做してよいことを仮定する. すなわち, 目的Seedの出現確率が消費数によらず常に$p$として良いものとする.

## 初期Seed厳選のモデリング

初期Seedから目的Seedまでに必要な消費数を$a$[**adv**]として, 閾値$U$[**adv**]に対して$a<U$を満たすならばその初期Seedを採択するという戦略を考える.

一回の試行にかかる時間(初期Seedを得るまでにかかる1回あたりの時間)を$C$[**s**], 乱数消費速度を$A$[**adv/s**]とする.

このとき, 試行回数$N$[**回**]で条件を満たす初期Seedを得られたのであれば, 最終的に目的Seedを得るまでにかかった時間$T$[**s**]は
\begin{align}
T = NC + \frac{a}{A} ~ \cdots(1)
\end{align}
と表される.

## "バラつき"の定量化
標本のバラつきを評価する指標には, "分散"と"標準偏差"の二種類が存在する. 分散は平均値からの偏差の2乗によって計算され, その平方根を取ることで標準偏差が得られる.

### "3シグマのルール"
平均を$\mu$, 標準偏差を$SD$とする正規分布において, 約$99.7\%$の
標本が$[\mu - 3SD, \mu + 3SD]$の中に含まれることが知られている.

また, 正規分布に従わない一般の分布の場合でも, 少なくとも$88.8\%$の標本は$[\mu - 3SD, \mu + 3SD]$の範囲内に入ることが知られている[^threesigma].
[^threesigma]:https://bigdata-analysis.hatenablog.com/entry/68_95_99-rule

すなわち, 目的のSeedを得るまでにかかる時間は概ね$E(T) + 3SD(T)$[**s**]で抑えられると言い換えられるので, これを最小化することを考える.

## 厳選時間の標準偏差の計算
$SD(T)$を計算するために, まずは分散$V(T)$を考える.

### 厳選時間の分散
$C, A$が定数と見なせることと試行回数$N$と必要な消費数$a$が独立に決まることに注意すると,

\begin{align}
V(T)
&= V(NC + \frac{a}{A})\\
&= V(N) \cdot C^2 + V(a) \cdot \frac{1}{A^2} ~ \cdots(2)
\end{align}

と変形される.

### 試行回数の分散
前回の考察にて, 初期Seed厳選が成功確率$p_{accept} = 1 - (1-p)^U$のベルヌーイ試行と見做せること,　ベルヌーイ試行が成功するまでの繰り返した際の試行回数$N$は幾何分布に従うことが分かっている.

幾何分布の分散は, ベルヌーイ試行の成功確率$p_{accept}$を用いて
\begin{align}
V(N) = \frac{1-p_{accept}}{p_{accept}^2}
\end{align}
として表される. $p, U$を用いて書き直すことで
\begin{align}
V(N)
&= \frac{1-(1 - (1-p)^U)}{(1 - (1-p)^U)^2}\\
&= \frac{(1-p)^U}{(1 - (1-p)^U)^2}
\end{align}
が導出される.

### 消費数の分散
前回の考察と同様に, $a$が$0$から$U-1$までの離散一様分布に従うと見做すことで,
\begin{align}
V(a) = \frac{U^2-1}{12}
\end{align}
が導出される.

### 総括
導出結果を纏めると,
\begin{align}
V(T) = \frac{C^2(1-p)^U}{(1 - (1-p)^U)^2} + \frac{U^2-1}{12A^2}
\end{align}
となる.

また, $|p|\ll 1$を仮定してマクローリン近似を行うことで,
\begin{align}
V(T)
&\approx \frac{C^2(1-pU)}{(1 - (1-pU))^2} + \frac{U^2-1}{12A^2}\\
&= \frac{C^2(1-pU)}{p^2U^2} + \frac{U^2-1}{12A^2}
\end{align}

と近似できる.従って, 目的Seedを得るまでにかかる時間の標準偏差
\begin{align}
SD(T) \approx \sqrt{\frac{C^2(1-pU)}{p^2U^2} + \frac{U^2-1}{12A^2}}
\end{align}
が計算できた.

## 目的関数の設計
導出された標準偏差と前回の期待値の近似の結果を用いて, 今回最小化を試みる目的関数$f(U)$を次のように定義する.
\begin{align}
f(U) &= \frac{C}{pU} + \frac{U}{2A} + 3\cdot \sqrt{\frac{C^2(1-pU)}{p^2U^2} + \frac{U^2-1}{12A^2}}
\end{align}

## 実験
パラメータを用意して$f(U)$を計算する. 今回はそれぞれ
- $C = 300$
- $A = 1.2$
- $p = \frac{1}{1000000}$

とした[^experiment].

[Wolfram Alphaのリンク](https://www.wolframalpha.com/input?i=f%28U%29+%3D+C%2F%28pU%29+%2B+U%2F%282A%29+%2B+3*sqrt%28C%5E2%281-pU%29%2F%28pU%29%5E2+%2B+%28U%5E2-1%29%2F%2812A%5E2%29%29+where+C%3D300%2C+A%3D1.2%2C+p%3D1%2F1000000&lang=ja)

[^experiment]:BDSPにおいて`g7TID`の1点狙いをする場合のモデル. 前回とは異なり, 自動化環境にて再測定した値を使用している.

### プロット
$f(U)$の計算式を [Wolfram Alphaに投入する](https://www.wolframalpha.com/input?i=f%28U%29+%3D+3+sqrt%28%2890000000000000000+%281+-+U%2F1000000%29%29%2FU%5E2+%2B+0.0578704+%28U%5E2+-+1%29%29+%2B+0.416667+U+%2B+300000000%2FU%2C+from+0+to+100000&lang=ja) と次のようになる.

### 最小値
前回と同様に, [Wolfram Alphaに最小値を計算させる.](https://ja.wolframalpha.com/input?i=%E6%A5%B5%E5%80%A4%E8%A8%88%E7%AE%97%E6%A9%9F&assumption=%22FSelect%22+-%3E+%7B%7B%22GlobalExtremaCalculator%22%7D%7D&assumption=%7B%22F%22%2C+%22GlobalExtremaCalculator%22%2C+%22curvefunction%22%7D+-%3E%223+sqrt%28%2890000000000000000+%281+-+U%2F1000000%29%29%2FU%5E2+%2B+0.0578704+%28U%5E2+-+1%29%29+%2B+0.416667+U+%2B+300000000%2FU%22&assumption=%7B%22F%22%2C+%22GlobalExtremaCalculator%22%2C+%22domain%22%7D+-%3E%220%3CU%22)

実験より, 目的関数$f(U)$の最小値を与える閾値$U$が求められた.

(以上, 常体終わり)

## 評価
前回の期待値の結果($U=26832$)と比較してみましょう.

[得られた目的関数に代入してみます.](https://www.wolframalpha.com/input?i=f%28U%29+%3D+C%2F%28pU%29+%2B+U%2F%282A%29+%2B+3*sqrt%28C%5E2%281-pU%29%2F%28pU%29%5E2+%2B+%28U%5E2-1%29%2F%2812A%5E2%29%29+where+C%3D300%2C+A%3D1.2%2C+p%3D1%2F1000000%2C+U%3D26832&lang=ja)

今回の最小値が, $f(U)\approx58708 ~(U\approx 32973)$ であったことを踏まえると, ざっくり$2000$[**s**]の短縮になったようです.

それなりに効果がありそうだと期待していたのですが, 思ってたほどではありませんでした.

とはいえ, 例えば[$f(100000)$の値をみる](https://www.wolframalpha.com/input?i=f%28U%29+%3D+C%2F%28pU%29+%2B+U%2F%282A%29+%2B+3*sqrt%28C%5E2%281-pU%29%2F%28pU%29%5E2+%2B+%28U%5E2-1%29%2F%2812A%5E2%29%29+where+C%3D300%2C+A%3D1.2%2C+p%3D1%2F1000000%2C+U%3D100000&lang=ja) と今回の最小値の二倍近く, 1日以上かかる可能性がある事が分かります.

この結果は, 最適戦略の重要性を示す一つの成果といえそうです.

## Appendix. 最小値を与える$\overline{U}$の推定
前回と同様に, $f(U)$の最小値を与えるような$\overline{U}$を考えたいのですが...
\begin{align}
f(U) &= \frac{C}{pU} + \frac{U}{2A} + 3\cdot \sqrt{\frac{C^2(1-pU)}{p^2U^2} + \frac{U^2-1}{12A^2}}
\end{align}

流石にこの式を解析的に解くのは難しそう...

一応, Excelの近似曲線機能と睨めっこしながらあれこれ試行錯誤したところ,
\begin{align}
\overline{U} \approx \sqrt{\frac{3CA}{p}}
\end{align}

が良さそうな近似を与えることが分かりました.

実際, 今回のパラメータ($C = 300$, $A = 1.2$, $p = \frac{1}{1000000}$)を代入すると,
\begin{align}
\overline{U}
&\approx \sqrt{\frac{3\cdot 300 \cdot 1.2}{\frac{1}{1000000}}}\\
&\approx 32863.4
\end{align}
となり, 真値に対する誤差が$\pm 0.5\%$ に収まっています.