あみだくじの横線は何本引くべきか

# 問題提起

一般的なあみだくじでは, 真下に行く確率が一番高くなり, 離れた箇所に行く確率が一番低い.
この性質は, 横線が少ないほど顕著であり, 多いと十分に無視できるようになるだろう.
では, 誤差が十分に無視できるような(十分に一様分布に近い)あみだくじを作るには何本横線が必要だろうか.

# 問題設定

あみだくじとして, 一番ベーシックな以下のパターンを考える.

* 縦線が計$N$本ある
* 横線が計$m$本ある
* 横線がどこを引くかは一様分布のランダムとする
* 横線は真横に引く

# 横線を一本だけ引いた場合
### ある場所nに一本引いた時
あみだくじの一番上の状態を$\vec{p}^{(0)}$とする.
仮に縦軸の右から$n$番目と$n+1$番目の間に横線を引いたとする.
そうしたときに$\vec{p}^{(0)}$の$n$番目の要素$p_n$と$n+1$番目の要素$p_{n+1}$を入れ替えるため, 
最終的な状態$\vec{p}^{(1)}$は
```math
\vec{p}^{(1)}=A_n p^{(0)}
\\
A_n=\left(\matrix{
1&& &\\\
&1 & &\\
&&\cdots&&\\
&&&0&1\\
&&&1&0\\
&&&&&\cdots\\
&&&&&&1&0\\
&&&&&&0&1
}\right )
```
と$n$,$n+1$番目だけ
```math
\matrix{0&1\\1&0}
```
にした単位行列$A_n$で表現できる.
この行列$A_n$は状態の遷移行列である.
($A_n$の$(i,j)$成分は状態$i$から状態$j$へ遷移する確率を返す)

### ランダムに横線を一本だけ引いたとき

横線の引き方は$1\to 2$から$(N-1) \to N$まで計$N-1$通り存在する.
そのため, ランダムな箇所に一本引いた際の遷移行列は, 各遷移行列$A_n$の期待値をとればいい.
```math
\vec{p}^{(1)}
=\frac{1}{N-1}\sum_n A_n \vec{p}^{(0)} \equiv A\vec{p}^{(0)}
\\
A=\frac{1}{N-1}\sum_n A_n=\left(\matrix{
1-\alpha & \alpha\\
\alpha & 1-2\alpha & \alpha\\
& \alpha & 1-2\alpha & \alpha\\
&&\cdots&\cdots&\cdots \\
&&&\alpha&1-2\alpha & \alpha\\
&&&&\alpha&1-\alpha
}\right)\\
\alpha=\frac{1}{N-1}
```

# ランダムに横線をm本引いた場合
横線を$m$本引いた場合のあみだくじも, 前節の考え方が適用できるため
```math
\vec{p}^{(m)}=A\vec{p}^{(m-1)}
```
このことから, 漸化的に計算ができ(マルコフ性)
```math
\vec{p}^{(m)}=A\vec{p}^{(m-1)}=A^2\vec{p}^{(m-2)}=\cdots =A^m\vec{p}^{(0)}
```
となる.
ここで,遷移行列$A$が$U^T\Lambda U$で対角化できるとすると
```math
\vec{p}^{(m)}=U^T\Lambda^mU\vec{p}^{(0)}
```
となる.
よって, $A$の固有値行列$\Lambda$を計算していく.

# 遷移行列Aの固有値計算
### 固有値の脚付き表記

$A$の要素$A_{ij}$は
```math
A_{ij}=\begin{cases}
1-\alpha & (i=j=1,N)
\\
1-2\alpha & (i=j\neq1,N)
\\
\alpha &(i=j\pm 1)
\\
0 &\mathrm{otherwise}
\end{cases}
```
固有値$\lambda$と固有ベクトル$\vec{u}$はその定義から
```math
0=(A-\lambda I)\vec{u}
```
そのため, 脚付き表記に書き直すと
```math
0=\sum_j (A_{ij}-\lambda \delta_{ij})u_j =\begin{cases}
\alpha u_{i+1}+(1-2\alpha -\lambda)u_i+\alpha u_{i-1} & (i\neq 1,N)
\\
(1-\alpha -\lambda)u_1 +\alpha u_2& (i=1)
\\
(1-\alpha -\lambda)u_N  +\alpha u_{N-1} & (i=N)
\end{cases}
```

### 隣接三項間漸化式を解く
計算簡略化のため, $u_0$と$u_{N+1}$分$u$を拡張する.
```math
\alpha u_{i+1}+(1-2\alpha -\lambda)u_i+\alpha u_{i-1}\ \ \ \ \ \ \ \ \ \ (1)
\\
u_0=u_1
\\
u_{N+1}=u_{N}
```

この３つのうち, 一番上は隣接三項間漸化式である.
そのため, 
```math
u_{i+1}+\gamma_\pm u_i =\gamma_\mp (u_i+\gamma_\pm u_{i-1})   \ \ \ \ \ \ \ \ \ \       (2)
```
という等比級数の形に式変形することで
```math
u_{i+1}+\gamma_\pm u_i =\gamma_\mp^{i} (1+\gamma_\pm )u_1 
```
と解くことができる.
このとき式1と式2の係数から$1=\gamma_+\gamma_-$が言えるため, $\gamma=\gamma_+=1/\gamma$とする.
また, 同様に係数の比較から
```math
\gamma+\gamma^{-1}=\frac{\lambda +2\alpha -1}{\alpha}
```
となる.

### λについて解く

一旦ここまでで, 使用できる数式をまとめる.
```math
u_{i+1}+\gamma^\pm u_i =\gamma^{\mp i-1} (1+\gamma^\pm )u_1 
\\
\gamma+\gamma^{-1}=\frac{\lambda +2\alpha -1}{\alpha}
\\
u_{N+1}=u_N
```
一番上の$(\pm,\mp)=(+,-)$の式と$(\pm,\mp)=(-,+)$の式の差から
```math
(\gamma -\gamma^{-1})u^i=[\gamma^{-i}(1+\gamma)-\gamma^{i}(1+\gamma^{-1})]u_1
```
$i=N$のとき$u_N=u_{N+1}$から
```math
\gamma^{-N}(1+\gamma)-\gamma^{N}(1+\gamma^{-1})=\gamma^{-N-1}(1+\gamma)-\gamma^{N+1}(1+\gamma^{-1})
```
これを解くと
```math
\gamma =\exp ( i\pi {n}/{N} )\ \ \ n=0,1,2,\cdots
```
となる.
このことから固有値$\lambda$は
$$
\lambda =1-2\alpha \left( 1-\frac{\gamma+\gamma^{-1}}{2} \right)=1-2\alpha\left\( 1-\cos \frac{\pi n}{2N} \right\)=1-4\alpha \sin^2\frac{\pi n}{2N}$$
$$(n=0,1,\cdots)$$

となる.
# あみだくじの横線は何本引くべきか
求めたい状態$\vec{p}^{(m)}$は固有値の$m$乗の線形結合で表現できる.
$$
    \vec{p}^{(m)}=\vec{q}^{(0)}+\lambda_1^m \vec{q}^{(1)}+\lambda_2^m \vec{q}^{(2)}+\cdots
$$
ここで, 第一項は$1-4\alpha \sin^2\frac{\pi \times 0}{2N}=1 $を用いた.
ここで, $\lambda_1,\lambda_2,\cdots $の絶対値は1未満のため, 
$$
\vec{p}^{(m)}\to \vec{q}_0 \ \ \ \ \ \ \ \ (m\to \infty) 
$$
である.
よって,  $\vec{q}_0$は一様分布と予想できる.
このことから, あみだくじに必要な横線の数$m$は,  $\lambda_1^m,\lambda_2^m,\cdots $が0とみなせる条件, すなわち二番目に大きい固有値$\lambda_1^m$が0とみなせる条件と言い換えられる.

仮に誤差が1%程度許容されると仮定すると, 引くべき横線の数$m$は
```math
\lambda_1^m <0.01
\\
m>\frac{\log 0.01}{\log \lambda}\simeq \frac{\log 0.01}{\pi^2} N^2(N-1)=0.467 N^2(N-1)
```
と見積もれる(詳しい計算はおまけ参照).
同様に誤差が5%程度許容されると仮定すると
```math
m>0.303 N^2(N-1)
```
# 数値計算での検証

実際に横線を$m$本引いた時の分布の一様分布のズレを求めた.

初期値は$p^{(0)}$をどこか一箇所が1のベクトルとして, 実際にあみだくじを繰り返すことで$p^{(m)}$を計算し
```math
RMSE=E[ (p^{(m)}-p_{一様})^2 ]^{1/2}
```
の指標を$m$ごとに調べた.

結果として, $0.467 N^2(N-1)$は誤差1%程度, $0.303 N^2(N-1)$は誤差5%程度に収束することが確認できた.
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2950973/7cd81a9d-ed73-6fb9-e6d9-7536530ac673.png)

# [おまけ]近似の計算
```math
\lambda_n=1-\frac{4}{N-1}\sin^2\frac{\pi n}{2N}=1-\frac{\pi^2 n^2}{N^2(N-1)}+\mathcal{O}(N^{-5})
\\
\log \lambda_n=-\sum_i \frac{1}{i}\left [1-\frac{\pi^2 n^2}{N^2(N-1)}+\mathcal{O}(N^{-5}) \right ]^i=-\frac{\pi^2 n^2}{N^2(N-1)}+\mathcal{O}(N^{-5})
```