しょぼんブログ

数学の色々とか様々とか

ドドスコードの出力回数の期待値

このプログラムにおいて, 配列{"ドド","スコ"}からランダムに要素を標準出力する回数の期待値を求めました.
また, そのあと「ドド」の出る確率を pとしたときの出力回数の期待値を求め,
「ドド」の出る確率を何にしたら期待値を最小にできるかを調べました.

今の状態を定める

もっと一般に, 「ドドスコスコスコ」が n回繰り返したものに一致するときに終了するものとします.
「ドド」と「スコ」をそれぞれ「 1個」として数えます.

「文字列 Sの末尾の i個が『ドドスコスコスコ× n』の開始 i個に一致して, どんな j>iに対してもSの末尾の j個が『ドドスコスコスコ× n』の開始 j個に一致しない」ことを, Sの状態は iであると言うことにします.
ここで, 任意の標準出力した文字列 Sに対して, その状態の値 iが一意に定まります.
（理由: 「 Sの末尾の i個が『ドドスコスコスコ× n』の開始 i個に一致する」という条件を満たす iの集合を Aとすると, Aは有限集合 \{0, 1, 2, \cdots, 4n\}の部分集合です. さらに 0 \in Aなので Aは空ではなく, よって最大元 iが存在します. このとき, Sは状態 iです.）

状態に対する期待値を定める

E_iを, 状態 iから終了するまでの出力回数の期待値とします. 文字列が異なっても, 上の状態が同じであれば, その状態から別の状態への遷移（後述）は共通します. なので, 状態が同じであるものをすべて同一視して, E_iを定めることができます.
状態 4nは最後の並びが「ドドスコスコスコ× n」となっているので, もう出力することはなくラブを注入することになります. よって E_{4n} = 0となります.
求めたいのは E_0で, これがドドスコードの出力回数の期待値になります.

E_nを求める

i \neq 4nとします.
iが 4で割り切れるとき, 標準出力 Sの末尾は「ドドスコスコスコ…ドドスコスコスコ」のようになっています. Sの末尾に「ドド」が加わった場合, 状態は i+1に進みます. 末尾に「スコ」が加わった場合, 残念ながら状態は 0に戻ってしまいます.
iが 4で割り切れないとき, 標準出力 Sの末尾は「ドドスコスコスコ…ドドスコ（スコが足りない）」のようになっています. Sの末尾に「ドド」が加わった場合, 残念ながら状態は 1に戻ってしまいます. 末尾に「スコ」が加わった場合, 状態は i+1に進みます.

各 iについて, E_iが満たす式は
E_i = \begin{cases} 0 & (i \geq 4n)\\1 + \frac{1}{2}(E_{i+1} + E_0) & (iは4で割り切れる)\\1 + \frac{1}{2}(E_{i+1} + E_1) & (iは4で割り切れない)\end{cases}
です. ここで, Sの末尾に「ドド」「スコ」が加わる確率は各 \frac{1}{2}で, さらに出力回数が +1されることに注意します.

たとえば, n=3（もとのツイートと同じ）のとき, 次のようになります.
E_{12} = 0
E_{11} = 1 + \frac{1}{2}(E_{12} + E_{1})
E_{10} = 1 + \frac{1}{2}(E_{11} + E_{1})
E_{9} = 1 + \frac{1}{2}(E_{10} + E_{1})
E_{8} = 1 + \frac{1}{2}(E_{9} + E_{0})
E_{7} = 1 + \frac{1}{2}(E_{8} + E_{1})
E_{6} = 1 + \frac{1}{2}(E_{7} + E_{1})
E_{5} = 1 + \frac{1}{2}(E_{6} + E_{1})
E_{4} = 1 + \frac{1}{2}(E_{5} + E_{0})
E_{3} = 1 + \frac{1}{2}(E_{4} + E_{1})
E_{2} = 1 + \frac{1}{2}(E_{3} + E_{1})
E_{1} = 1 + \frac{1}{2}(E_{2} + E_{1})
E_{0} = 1 + \frac{1}{2}(E_{1} + E_{0})
これは 4n+1元の連立一次方程式になります. 行列を考えれば, 掃き出し法によって解くこともできますが, E_0, E_1の値がわかっていれば, 後ろから順番に求めていくことで E_2から E_{4n}まですべて求まります. なので, E_0, E_1を文字で置いて仮定することを考えます.

E_0, E_1を仮定する

E_0 = \alpha, E_1 = \betaとおきます. ここで, E_i = N_i + A_i \alpha + B_i \betaとして数列 N_i, A_i, B_iをおきます.

これをふまえて, E_iは
E_i = \begin{cases} 0 & (i \geq 4n)\\1 + \frac{1}{2}(N_{i+1} + A_{i+1} \alpha + B_{i+1} \beta + \alpha) & (iは4で割り切れる)\\1 + \frac{1}{2}(N_{i+1} + A_{i+1} \alpha + B_{i+1} \beta + \beta) & (iは4で割り切れない)\end{cases}

よって, N_iについては
N_i = \begin{cases} 0 & (i \geq 4n)\\1 + \frac{1}{2}N_{i+1} & (他)\end{cases}
A_iについては
A_i = \begin{cases} 0 & (i \geq 4n)\\\frac{1}{2}A_{i+1} + \frac{1}{2} & (iは4で割り切れる)\\\frac{1}{2}A_{i+1} & (iは4で割り切れない)\end{cases}
B_iについては
B_i = \begin{cases} 0 & (i \geq 4n)\\\frac{1}{2}B_{i+1} & (iは4で割り切れる)\\\frac{1}{2}B_{i+1} + \frac{1}{2} & (iは4で割り切れない)\end{cases}
となります.

これは, 簡単な漸化式になっています！ N_2, A_2, B_2は, 簡単な式で表すことができます（後述）. N_2, A_2, B_2を具体的に求めるまえに, それらで \alpha, \betaをどうやって解くかを書きます.

N_2, A_2, B_2が分かったとき

E_{1} = 1 + \frac{1}{2}(E_{2} + E_{1})
E_{0} = 1 + \frac{1}{2}(E_{1} + E_{0})
が成立するので,
\beta = 1 + \frac{1}{2}(N_2 + A_2 \alpha + B_2 \beta + \beta)…①
\alpha = 1 + \frac{1}{2}(\beta + \alpha)…②
となります.

②を変形すると, \alpha = 2 + \beta…②'となります.
①を変形すると, \beta = 1 + \frac{N_2}{2} + \frac{A_2}{2} \alpha + (\frac{B_2}{2} + \frac{1}{2}) \beta
2倍して②'を使うと 2 \beta = 2 + N_2 + 2A_2 \alpha + (A_2 + B_2 + 1) \betaとなり, \betaについて解くと \displaystyle \beta = \frac{2+N_2 + 2A_2}{1-A_2-B_2}がわかります.

なので, 求める期待値は \displaystyle E_0 = \alpha = \beta + 2 = \frac{2+N_2 + 2A_2}{1-A_2-B_2} + 2です.

1-A_2-B_2, N_2を求める

先ほどの

A_i = \begin{cases} 0 & (i \geq 4n)\\\frac{1}{2}A_{i+1} + \frac{1}{2} & (iは4で割り切れる)\\\frac{1}{2}A_{i+1} & (iは4で割り切れない)\end{cases}
B_i = \begin{cases} 0 & (i \geq 4n)\\\frac{1}{2}B_{i+1} & (iは4で割り切れる)\\\frac{1}{2}B_{i+1} + \frac{1}{2} & (iは4で割り切れない)\end{cases}

式より, A_i + B_iをまとめて考えると
A_i+B_i = \begin{cases} 0 & (i \geq 4n)\\\frac{1}{2}(A_{i+1}+B_{i+1}) + \frac{1}{2} & (他)\end{cases}
となります.
これは簡単に解くことができます. 特性方程式を考えるなどして, A_i+B_i - 1 = \frac{1}{2}(A_{i+1}+B_{i+1} - 1)が得られます. 数列(A_i+B_i-1)は等比数列なので, A_2+B_2-1 = (-1) \times (\frac{1}{2})^{4n-2}を得ます. よって, 1-A_2-B_2 = \frac{1}{2^{4n-2}}です.

N_2を求めてみましょう.
N_i = \begin{cases} 0 & (i \geq 4n)\\1 + \frac{1}{2}N_{i+1} & (他))\end{cases}
です. 同様に N_i-2 = \frac{1}{2}(N_{i+1} - 2)より N_2-2 = (-2) \times (\frac{1}{2})^{4n-2} よって N_2 = 2 - \frac{2}{2^{4n-2}} となります.

A_2を求める

A_iは iが 4の倍数のときだけに新たに \frac{1}{2}を生み出し, それ以外は後ろの項を \frac{1}{2}に減衰させているというイメージを持つことができます. よって, A_{4i}に注目すると, \displaystyle A_{4i} = \frac{1}{2} A_{4i+1} + \frac{1}{2} = \frac{1}{4} A_{4i+2} + \frac{1}{2} = \frac{1}{8} A_{4i+3} + \frac{1}{2} = \frac{A_{4(i+1)}}{16} + \frac{1}{2}です.
ここにくると, 同様に \displaystyle A_{4i} - \frac{8}{15} = \frac{1}{16} \left( A_{4(i+1)} - \frac{8}{15}\right) より \displaystyle A_4 = \frac{8}{15}\left(1-\frac{1}{16^{n-1}}\right) = \frac{2^{4n-4} - 1}{15\times 2^{4n-7}} が従います.

A_2 = \frac{1}{2}A_3 = \frac{1}{4}A_4なので, \displaystyle A_2 = \frac{2^{4n-4} - 1}{15\times 2^{4n-5}} です.

ついに E_0 が求められるぞ

材料は揃いました. E_0 = \alphaを求めましょう. さきほど, \displaystyle E_0 = \frac{2+N_2 + 2A_2}{1-A_2-B_2} + 2 を求めました. これに,
\displaystyle N_2 = 2 - \frac{2}{2^{4n-2}}
\displaystyle A_2 = \frac{2^{4n-4} - 1}{15\times 2^{4n-5}}
\displaystyle 1-A_2-B_2 = \frac{1}{2^{4n-2}}
を代入すると,

\displaystyle E_0 = 2^{4n-2} \times \left(2 + 2 - \frac{2}{2^{4n-2}} + \frac{2^{4n-4} - 1}{15\times 2^{4n-5}}\right) + 2\\ \displaystyle = 2^{4n-2} \times \frac{4 \times 2^{4n-2} - 2 + 2\times 8 (2^{4n-4} - 1)}{15 \times 2^{4n-2}} + 2\\ \displaystyle = 4 \times 2^{4n-2} - 2 + 2\times \frac{16 (2^{4n-4} - 1)}{15} + 2\\ \displaystyle = 16 \times 2^{4n-4} + \frac{16}{15} \times (2^{4n-4} - 1)\\ \displaystyle = \frac{16\times(15+1) 2^{4n-4} - 16}{15}\\ \displaystyle = \frac{16(2^{4n} - 1)}{15}
となり, \displaystyle E_0 = \frac{16(2^{4n} - 1)}{15}が導かれました. 思ったより簡単な式になりました.

これに n =3を代入したものが, 目的の期待値になります. その期待値は, \displaystyle \frac{16(4096 - 1)}{15} = 4368 です.
なので, プログラムを動かすと, 平均で 4368回ものドドスコを覚悟しなければなりません.
2^{4n}-1 = 16^n-1 は合同式などを考えることによって, 15で割り切れることがわかります. よって, nがどんな正整数でもドドスコードの出力回数の期待値は整数になります.

以下は n=1, \cdots, 10の期待値です.
n = 1 \to E_0 = 16\\ n = 2 \to E_0 = 272\\ n = 3 \to E_0 = 4368\\ n = 4 \to E_0 = 69904\\ n = 5 \to E_0 = 1118480\\ n = 6 \to E_0 = 17895696\\ n = 7 \to E_0 = 286331152\\ n = 8 \to E_0 = 4581298448\\ n = 9 \to E_0 = 73300775184\\ n = 10 \to E_0 = 1172812402960\\
出力回数の分布を調べてないので詳しいことは言えませんが, n\geq 10でドドスコードを実行するのは危険そうです.

偏りを持ったドドスコ

さっきまでは「ドド」「スコ」は均等に出現するときを考えていました. 今度は, 「ドド」と「スコ」で偏りがある場合を考えてみます.「ドド」を出す確率を p ~ (0\lt p\lt 1), 「スコ」を出す確率を 1-pとします.
同様に考えると,
E_iは
E_i = \begin{cases} 0 & (i \geq 4n)\\1 + pE_{i+1} + (1-p)E_0 & (iは4で割り切れる)\\ 1 + pE_1 + (1-p)E_{i+1} & (iは4で割り切れない)\end{cases}
よって,
N_i = \begin{cases} 0 & (i \geq 4n)\\1 + pN_{i+1}& (iは4で割り切れる)\\ 1 + (1-p)N_{i+1} & (iは4で割り切れない)\end{cases}
A_i = \begin{cases} 0 & (i \geq 4n)\\(1-p) + pA_{i+1} & (iは4で割り切れる)\\(1-p)A_{i+1} & (iは4で割り切れない)\end{cases}
B_i = \begin{cases} 0 & (i \geq 4n)\\pB_{i+1} & (iは4で割り切れる)\\(1-p)B_{i+1} + p & (iは4で割り切れない)\end{cases}
となります.

とくに,
E_1 = 1 + pE_1 + (1-p)E_2
E_0 = 1 + pE_1 + (1-p)E_0
であるので, 同様に連立方程式を考えると
\displaystyle \alpha = \frac{1}{p} + \beta より \displaystyle \beta = 1+ (1-p)N_2 + \frac{1-p}{p}A_2 + (p+(1-p)(A_2+B_2)) \beta
すなわち,
\displaystyle \alpha = \frac{1}{p} + \frac{1}{1-A_2-B_2} \left(\frac{1}{1-p} + N_2 + \frac{1}{p} A_2\right)
となります.

さきほど A_2を求めるときに使ったテクニックを使います.

N_iについて N_{4i} = 2-(1-p)^3 + (1-p)^3 p N_{4i+4} より \displaystyle N_4 = \frac{2-(1-p)^3}{1-(1-p)^3p} \left(1 - \{(1-p)^3p\}^{n-1}\right) です.
よって, \displaystyle N_2 = 2-p + (1-p)^2 \frac{2-(1-p)^3}{1-p(1-p)^3} \left(1 - \{p(1-p)^3\}^{n-1}\right) です.

A_iについて A_{4i} = (1-p) + p(1-p)^3A_{4i+4} より \displaystyle A_4 = \frac{1-p}{1-p(1-p)^3} (1-\{p(1-p)^3\}^{n-1}) です.
よって, \displaystyle A_2 = \frac{(1-p)^3}{1-p(1-p)^3} (1-\{p(1-p)^3\}^{n-1})です.

A_i + B_iについて A_{4i}+B_{4i} = 1-p(1-p)^3 + p(1-p)^3 (A_{4i+4} + B_{4i+4}) より \displaystyle A_4 + B_4 = 1-(p(1-p)^3)^{n-1} です.
よって, \displaystyle A_2 + B_2 = 2p-p^2 + (1-p)^2 (1-\{p(1-p)^3\}^{n-1}) です.

さあ, \alphaを求めていきます.
まず, 1-A_2-B_2 = 1-2p+p^2 -(1-p)^2 (1-\{p(1-p)^3\}^{n-1})\\ = (1-p)^2 \{1 - (1-\{p(1-p)^3\}^{n-1})\} = (1-p)^2 \{p(1-p)^3\}^{n-1}

です.よって,
\displaystyle E_0 = \alpha = \frac{1}{p} + \frac{1}{1-A_2-B_2} \left(\frac{1}{1-p} + N_2 + \frac{1}{p} A_2\right)\\ \displaystyle = \frac{1}{p} + \frac{\frac{1}{1-p}+2-p + (1-p)^2 \frac{2-(1-p)^3}{1-p(1-p)^3} \left(1 - \{p(1-p)^3\}^{n-1}\right)}{(1-p)^2 \{p(1-p)^3\}^{n-1}}\\ \displaystyle + \frac{\frac{(1-p)^3}{p(1-p(1-p)^3)} (1-\{p(1-p)^3\}^{n-1})}{(1-p)^2 \{p(1-p)^3\}^{n-1}}\\ \displaystyle = \frac{1}{p} + \frac{p^2-3p+3}{(1-p)^3\{p(1-p)^3\}^{n-1}} + \frac{(1+p-p(1-p)^3)(1-\{p(1-p)^3\}^{n-1})}{p(1-p(1-p)^3)\{p(1-p)^3\}^{n-1}}

です. おぞましいです.

期待値が最小となる pの値

これ以上は厳しそうなので, Desmosさんに助けてもらいます.
https://www.desmos.com/calculator/fydiaq65et
↑は各 pに対する期待値の対数をとったものを表示させています. こんなグラフも対応できるDesmosさんはすごいです.

ということで, n = 3のときの期待値の最小値は, p=0.25のとき, すなわち「ドド」が25%, 「スコ」が75%のときの951.751らしいです. 最速でラブを注入したいときはこの比率に設定しましょう.
ちなみに, 実験の結果によると, どんな nでも p=0.25のとき最小になりそうです.（証明はできていません） 直感的には0.25周辺になる感じはありますが, 0.25ピッタリというとちょっとピンと来ません.

以下は n=1, \cdots, 10の p=0.25での期待値です.（小数点第2位まで）
n = 1 \to E_0 = 9.48\\ n = 2 \to E_0 = 99.38\\ n = 3 \to E_0 = 951.75\\ n = 4 \to E_0 = 9033.49\\ n = 5 \to E_0 = 85660.35\\ n = 6 \to E_0 = 812196.46\\ n = 7 \to E_0 = 7700835.17\\ n = 8 \to E_0 = 73015335.57\\ n = 9 \to E_0 = 692293561.60\\ n = 10 \to E_0 = 6563968593.52
こっちは n=10でもいけそうです.

p=0.25のときの期待値を e(0.25), p=0.50のときの期待値を e(0.50)とするとき, e(0.50)/e(0.25)の値は以下のようになります.（小数点第2位まで）
n = 1 \to e(0.50)/e(0.25) = 1.69\\ n = 2 \to e(0.50)/e(0.25) = 2.74\\ n = 3 \to e(0.50)/e(0.25) = 4.59\\ n = 4 \to e(0.50)/e(0.25) = 7.74\\ n = 5 \to e(0.50)/e(0.25) = 13.06\\ n = 6 \to e(0.50)/e(0.25) = 22.03\\ n = 7 \to e(0.50)/e(0.25) = 37.18\\ n = 8 \to e(0.50)/e(0.25) = 62.74\\ n = 9 \to e(0.50)/e(0.25) = 105.88\\ n = 10 \to e(0.50)/e(0.25) = 178.67
nが大きいほど p=0.25の方がいいことがわかります.

n=10で試してみました. 標準出力をいっぱいするのは厳しいので, 擬似的なものになりますが…

import random
from collections import deque
import datetime

print(datetime.datetime.now())

n = 10
a = deque()
b = deque()
for i in range(4*n):
	if i % 4 == 0:
		b.append("ドド")
	else:
		b.append("スコ")
	a.append("")

i = 0
while True:
	a.popleft()
	if random.randrange(4) == 0:
		a.append("ドド")
	else:
		a.append("スコ")
	i += 1
	if a == b:
		print(i)
		break
	if i % 100000000 == 0:
		print(datetime.datetime.now(), i)

a.append("ラブ注入♡")
print(*a,sep="")

なんと期待値の1/20で帰らせてくれました. 今回はやさしかったです.

感想

ぼくは大学の必修科目の課題がまだ残っているにもかかわらず, ドドスコの出力回数の期待値を求めるのに心血を注ぎました　は？
でも数字をいじるのはたのしかったです