さいころの目の和がkの倍数になる確率

　どうも, しょぼんです. この記事では$m$面さいころ（特に$m=6$）を$n$回投げて出た目の総和が$k$の倍数になる確率について書きます. 学部1年生で習う線形代数の知識を使います.

この記事で扱う結果

1. 任意の$m$での一般項を導き, 特に$m=6$で$k=1,2,3,4,5,6,7,8$の場合の一般項を求める.
2. $m\ge 2$なら, どのような$k$でも$n\to \mathrm{\infty }$のとき確率は$1/k$に収束する.

漸化式を立てる

　$n$回投げて$k$で割った余りが$l$になる確率を$p[n,l]$とします. このとき, $n+1$回投げて$k$で割った余りが$l$になる確率は, $a\mathrm{\%}b$を$a$を$b$で割った余りとして

• $k$で割って余りが$(l-1)\mathrm{\%}k$のときから$1/m$の確率で$1$が出る
• $k$で割って余りが$(l-2)\mathrm{\%}k$のときから$1/m$の確率で$2$が出る
• $⋮$
• $k$で割って余りが$(l-m)\mathrm{\%}k$のときから$1/m$の確率で$m$が出る
  場合だけなので, 次の漸化式が立てられます.
  $p[n+1,l]=\frac{1}{m}(p[n,(l-1)\mathrm{\%}k]+p[n,(l-2)\mathrm{\%}k]+\cdots +p[n,(l-m)\mathrm{\%}k])$
  ここで, $p[0,0]=1,p[0,1]=\cdots =p[0,k-1]=0$とします. 理由は何もたさない場合, その和は$0$と考えるからです.
  　たとえば$6$面さいころを投げて$4$の倍数になる確率を考える場合,
  $$\begin{gathered} p[n+1,0]=\frac{1}{6}(p[n,0]+p[n,1]+2p[n,2]+2p[n,3]) \\ p[n+1,1]=\frac{1}{6}(2p[n,0]+p[n,1]+p[n,2]+2p[n,3]) \\ p[n+1,2]=\frac{1}{6}(2p[n,0]+2p[n,1]+p[n,2]+p[n,3]) \\ p[n+1,3]=\frac{1}{6}(p[n,0]+2p[n,1]+2p[n,2]+p[n,3]) \end{gathered}$$
  となり, 行列で表すと,
  ${\displaystyle \left(\begin{array}{c}p[n+1,0]\\ p[n+1,1]\\ p[n+1,2]\\ p[n+1,3]\end{array}\right)=\frac{1}{6}\left(\begin{array}{cccc}1& 1& 2& 2\\ 2& 1& 1& 2\\ 2& 2& 1& 1\\ 1& 2& 2& 1\end{array}\right)\left(\begin{array}{c}p[n,0]\\ p[n,1]\\ p[n,2]\\ p[n,3]\end{array}\right)}$
  となります.

特別な形の$k$について

　$k$の値によって, 漸化式がかんたんに解ける場合があります.

$k$が$m$の約数である場合

　たとえば$6$面サイコロを$n$回振って$3$の倍数になる場合などです. この場合, $m/k$は整数なので$(l-1)\mathrm{\%}k,(l-2)\mathrm{\%}k,\cdots ,(l-m)\mathrm{\%}k$に対し$0,1,\cdots ,k-1$はどれも均等に$m/k$回出てきます. よって

さて、$p[n,0]+p[n,1]+\cdots +p[n,k-1]=1$なので

となります. 上の漸化式は$n\ge 0$に成り立つので, $p[n,0]$は$n\ge 1$のとき$p[n,0]=\frac{1}{k}$になります.

$m=6$の場合の結果

　$6$の正の約数は$1,2,3,6$であるので, 上の議論が使えて, $n\ge 1$とすると

$k=1$のとき$p[n,0]=1$
$k=2$のとき$p[n,0]=\frac{1}{2}$
$k=3$のとき$p[n,0]=\frac{1}{3}$
$k=6$のとき$p[n,0]=\frac{1}{6}$

ということがわかりました.

$k=m+1$の場合

　たとえば$6$面サイコロを$n$回振って$7$の倍数になる場合などです. この場合, $(l-1)\mathrm{\%}k,(l-2)\mathrm{\%}k,\cdots ,(l-m)\mathrm{\%}k$は$0,1,\cdots ,k-1$のうち$l$以外が$1$回ずつ出てきます. よって

となって, $l$以外の状態に依存しません. これは隣接2項間漸化式となり, 解くことができます. 高校数学でよくある形です.
　$\alpha =\frac{1}{m}-\frac{1}{m}\alpha$を解くと$\alpha =\frac{1}{m+1}$となるので,

$p[n,l]-\frac{1}{m+1}=q[n]$とすると, $q[n+1]=-\frac{1}{m}q[n]$であり等比数列になるので, $q[n]=(-\frac{1}{m}{)}^{n}q[0]$とわかります. $l=0$のとき$q[0]=\frac{m}{m+1}$なので

がわかりました. $l\ne 0$のときは$q[0]=-\frac{1}{m+1}$となるので$p[n,l]=\frac{1}{m+1}-\frac{1}{m+1}(-\frac{1}{m}{)}^n$となります.

$m=6$の場合の結果

　$k=7$に対して上の議論が使えて,${\displaystyle p[n,0]=\frac{1}{7}\{1-{(-\frac{1}{6})}^{n-1}\}}$となります.

一般の$k$について

　いま, $k$が特別な形のときを解きましたが, 一般の$k$に対してはどうなるでしょうか？たとえば$m=6$のとき$k=4,5,8$などです.

行列累乗で表す

　$p[n+1,l]$は$p[n,0],p[n,1],\cdots ,p[n,k-1]$の一次結合で表されるので

を満たす$k$次正方行列$A$があります. これを繰り返し使うと,$n\ge 1$のとき

となります. もし$A$が行列$P$によって対角行列$D$に対角化される, すなわち$P^{-1}AP=D$のとき, $A^n=P{D}^n{P}^{-1}$となります. $A$の固有値を${\alpha }_0,{\alpha }_1\cdots ,{\alpha }_{k-1}$として

とするとき,

となるので, $P,P^{-1},D$が求まれば行列の積を$2$回行うことで簡単に一般項が計算できます.

$A$は巡回行列になる

　巡回行列とは,

と表される行列のことです. $A$が巡回行列となることを示しましょう. まず,

と表されるとします.

より,

ということで

となり$A$は巡回行列です. たとえば$m=6,k=4$のとき

となることは先ほど求めましたが, $c_0=\frac{1}{6},c_1=\frac{1}{6},c_2=\frac{2}{6},c_3=\frac{2}{6}$ として

が確かめられます.

巡回行列の固有値と固有ベクトル

　巡回行列

について, $\omega =\exp \left(\frac{2\pi i}{n}\right)=\cos \left(\frac{2\pi }{n}\right)+i\sin \left(\frac{2\pi }{n}\right)$, $f(x)=c_0+c_{1}x+\cdots +c_{n-1}{x}^{n-1}=\sum _{k=0}^{n-1}{c}_k{x}^k$として$n$つの複素数$f(1),f(\omega ),\cdots ,f({\omega }^{n-1})$が固有値になり, それぞれに対応する固有ベクトルは$c_k$の値に関係なく

となることが知られています. さらにそれらは標準内積で正規直交基底になります. これは直接示せます. $0\le k<n$を整数として,

の第$m$成分を比較すると, 左辺は$\frac{1}{\sqrt{n}}\sum _{k=0}^{n-1}{c}_{(-m+k)\mathrm{\%}n}{\omega }^k$, 右辺は$\frac{1}{\sqrt{n}}\sum _{k=0}^{n-1}{c}_k{\omega }^{k+m}$です. $\omega$は${\omega }^n$を満たすことに注意すれば

（ただし$m=0$の場合は第2項は省略）となって一致します. さらに固有値$f({\omega }^k),f({\omega }^m)$に対する固有ベクトルの標準内積の第$l$項について, $\frac{1}{n}(1+\overline{{\omega }^k}{\omega }^m+\cdots +\overline{{\omega }^{k(n-1)}}{\omega }^{m(n-1)})=\frac{1}{n}(1+{\omega }^{m-k}+\cdots +{\omega }^{(m-k)(n-1)})$となるので, $m=k$のとき内積は$1$, そうでないなら内積は$0$となります. よって固有ベクトル

は正規直交基底をなします.
　以上の性質から, $n$次の巡回行列なら

として$P$は正規直交基底であることから$P$はユニタリ行列（${}^t\bar{P}P$が単位行列$E$になる）で

を満たします. さらに,

と対角化できます. ちなみに$P$はDFT行列という名前があり, 離散フーリエ変換と関係があるそうです（詳しくないので今は書けません…）

一般の$k$について$p[n,0]$を求める

　いよいよ$m$面さいころを$n$回投げて出た目の総和が$k$の倍数になる確率を求めます. $k$次の巡回行列$A$によって

と表されるので, 上の節で見た通り, $D=P^{-1}AP$で対角化されます. $A^n=P{D}^n{P}^{-1}$より,

となります. 一番うしろのベクトルが第$1$成分以外$0$であることに注意すれば, $P{D}^n{P}^{-1}$の$1$列目だけを見ればよいです.

で,

より

となり, ${\displaystyle p[n,0]=\frac{1}{k}(f(1{)}^n+f(\omega {)}^n+\cdots +f({\omega }^{k-1}{)}^n)}$ということがわかりました. さらに$c_0+c_1+\cdots +c_{k-1}=1$を考慮すれば, $f(1)=1$ということがわかるので, ${\displaystyle p[n,0]=\frac{1}{k}(1+f(\omega {)}^n+\cdots +f({\omega }^{k-1}{)}^n)}$ となります.

$m=6,k=4$のとき

となるので, $f(x)=\frac{1}{6}(1+x+2x^2+2x^3)$ となります. $\omega =\exp \left(\frac{i\pi }{2}\right)=\cos \left(\frac{\pi }{2}\right)+i\sin \left(\frac{\pi }{2}\right)=i$より
$f(\omega )=\frac{1}{6}(1+i+2i^2+2i^3)=\frac{1}{6}(-1-i)$
$f({\omega }^2)=\frac{1}{6}(1+i^2+2i^4+2i^6)=0$
$f({\omega }^3)=\frac{1}{6}(1+i^3+2i^6+2i^9)=\frac{1}{6}(-1+i)$
となります. $f(\omega )=\frac{\sqrt{2}}{6}(\cos (-\frac{3\pi }{4})+i\sin (-\frac{3\pi }{4})),f({\omega }^3)=\frac{\sqrt{2}}{6}(\cos \left(\frac{3\pi }{4}\right)+i\sin \left(\frac{3\pi }{4}\right))$ より,

となって, ${\displaystyle p[n,0]=\frac{1}{4}(1+f(\omega {)}^n+f({\omega }^2{)}^n+f({\omega }^3{)}^n)}$だったので,

が導けました！いくつか項を計算すると,
$p[1,0]=\frac{1}{6^1}=\frac{1}{6}=0.1666\cdots$
$p[2,0]=\frac{9}{6^2}=\frac{1}{4}=0.25$
$p[3,0]=\frac{55}{6^3}=\frac{55}{216}=0.2546\cdots$
$p[4,0]=\frac{322}{6^4}=\frac{161}{648}=0.2484\cdots$
$p[5,0]=\frac{1946}{6^5}=\frac{973}{3888}=0.2502\cdots$
となります.

今回は$p[n,1]$も計算してみます. ${\displaystyle p[n,1]=\frac{1}{4}(1+\omega f(\omega {)}^n+{\omega }^{2}f({\omega }^2{)}^n+{\omega }^{3}f({\omega }^3{)}^n)}$ですが

となり

です. 同様にして

となります. 綺麗ですね！

$m=6,k=5$のとき

より, $f(x)=\frac{1}{6}(1+x+x^2+x^3+2x^4)$ となります. $\omega =\exp \left(\frac{2i\pi }{5}\right)=\cos \left(\frac{2\pi }{5}\right)+i\sin \left(\frac{2\pi }{5}\right)$であるので,

に注意して
$f(\omega )=\frac{1}{6}{\omega }^4=\frac{1}{6}{\omega }^{-1}=\frac{1}{6}(\cos \left(\frac{-2\pi }{5}\right)+i\sin \left(\frac{-2\pi }{5}\right))$
$f({\omega }^2)=\frac{1}{6}({\omega }^2{)}^4=\frac{1}{6}{\omega }^3=\frac{1}{6}{\omega }^{-2}=\frac{1}{6}(\cos \left(\frac{-4\pi }{5}\right)+i\sin \left(\frac{-4\pi }{5}\right))$
$f({\omega }^3)=\frac{1}{6}({\omega }^3{)}^4=\frac{1}{6}{\omega }^2=\frac{1}{6}(\cos \left(\frac{4\pi }{5}\right)+i\sin \left(\frac{4\pi }{5}\right))$
$f({\omega }^4)=\frac{1}{6}({\omega }^4{)}^4=\frac{1}{6}\omega =\frac{1}{6}(\cos \left(\frac{2\pi }{5}\right)+i\sin \left(\frac{2\pi }{5}\right))$
となるので,

となって, ${\displaystyle p[n,0]=\frac{1}{5}(1+f(\omega {)}^n+f({\omega }^2{)}^n+f({\omega }^3{)}^n+f({\omega }^4{)}^n)}$だったので,

が導けました！いくつか項を計算すると,
$p[1,0]=\frac{1}{6^1}=\frac{1}{6}=0.1666\cdots$
$p[2,0]=\frac{7}{6^2}=\frac{7}{36}=0.1944\cdots$
$p[3,0]=\frac{43}{6^3}=\frac{43}{216}=0.1990\cdots$
$p[4,0]=\frac{259}{6^4}=\frac{259}{1296}=0.1998\cdots$
$p[5,0]=\frac{1556}{6^5}=\frac{389}{1944}=0.2001\cdots$
となります. 同様に$p[n,1]$から$p[n,4]$も計算してみると