コイン投げから一様分布（と正規分布）を生成する方法

表題の通り
ファイナンスのための確率解析Ⅱ（S.E.シュリーブ）の第一章に書いてあった一様分布の生成方法を検証してみる。

具体的な方法

表が出る確率が1/2であるコインを考える。（要はp = 1/2のベルヌーイ試行）
このコインをn回投げるとする。n = 1,2,...に対して
\[
Y_n(\omega) = \begin{cases}
1 & \omega_n = H \\
0 & \omega_n = T
\end{cases}
\]
そして、
\[
X = \sum^{\infty}_{n=1} \frac{Y_n}{2^n}
\]
と定義する。

この時、Xは[0,1]上に一様に分布する。

検証

MATLABを用いて一様分布に従っているか調べる。残念ながら私は高々有限回しかコインを投げられないので*1、とりあえずn = 5,10,100と値を変えて調べることにする。また、サンプルサイズは1000とする。

結果

階級の幅を0.1としたヒストグラムは以下の通りである。

また、縦軸に経験分布、横軸を理想的な一様分布としたqqプロット*2は以下の通りである。

qqプロット

プロットされたデータを見ると、n=5でも十分に一様分布しているように見える。（一方でnを増やすと、より一様分布するとは限らない？）

また、定性的な結果を得るために上のヒストグラムを作成した際に離散化した乱数値に対してカイ二乗適合度検定*3を行う。自由度は９である。
有意水準5％において、すべてのnに対して帰無仮説は棄却されなかった。従って、離散化した乱数値は離散一様分布に従っていないとは言えない*4ということが分かった。*5

更なる検証

5回コインを投げればそれなりの一様乱数を生成できることが分かったので、これを用いて標準正規乱数を作成してみる。ここでは、ボックス＝ミュラー法*6を用いることにする。
尚、ボックス＝ミュラー法では一様乱数が二つ必要なので、乱数は二倍必要とされる。（n=5であれば10回、n=10であれば20回コインを投げる必要がある。）

結果は以下の通りである。

ヒストグラムを見る限りでは正規分布に従っているようだが、qqプロットによると、裾の方が標準正規分布から乖離していることがわかる。

また、生成された乱数に対してコルモゴロフースルミノフ検定*7を行ったところ、帰無仮説は棄却された。従って、標準正規乱数であるとは言えないことが分かった。

nの値を増やしてもう一度検証をしてみる。とりあえずn=10とする。

今度はqqプロットの結果も含め正規分布に従っているように見える。実際、コルモゴロフ-スルミノフ検定の結果は棄却されなかった。

まとめ

コイン投げの結果から一様乱数、更には標準正規乱数を生成してみようと試みた。結果としては、単に一様分布を生成したいだけであればn=5でもある程度十分な精度の乱数が得られる事が分かった。
ただ、そこから更に標準正規乱数を生成しようとするとn=5では不十分であり、n=10、つまり20回コインを投げる必要がある事が分かった。(20回もコイン投げるんだったら二項分布から正規分布への近似法使った方がよくね？？）

もし無人島に裸一貫で飛ばされたとしても、（歪んでいない）コインが一枚あれば一様分布、正規分布、さらには逆関数法を使えば他にも様々な分布に従う乱数が作れると思うと夢が膨らみますね！！（強引な結論）

使用したコード（matlab）

%コイン投げから一様分布を生成し、プロットや検定を行うコード
n_list = [5, 10, 100] %ループ回数のリスト

for i = 0:2
    
    n = n_list(1, i + 1) %ループ回数を設定
    
    sample_size = 1000;
    
    %乱数の配列を作る（sample_size × 1)
    ransu = double.empty(sample_size,0); %初期配列
    for j = 1:sample_size
        ransu(j,1) = itiyo(n);
    end
    
　%ヒストグラムを作成する
    figure
    histogram(ransu)
    xlabel("乱数の値")
    ylabel("度数")
    title(['乱数の分布 n = ', num2str(n)], 'FontSize',12)
    
    %乱数が一様分布に従っているか、カイ二乗適合度検定を行う。(自由度は9とする。）
    %%カイ二乗統計量を求める
    toukeiryo = kai2_kentei(sample_size, ransu)
    
    %%上側確率を求める。
    upper_prob = 1 - chi2cdf(toukeiryo,9)
    
    %%結果をansに入れる。
    if upper_prob >= 0.05
        ans = "一様分布していないとはいえない" %棄却できない
    else
        ans = "一様分布していない" %棄却できる
    end
    
    %qqプロットを行う
    pd = makedist('Uniform');
    
    figure;
    qqplot(ransu,pd);
    title(['標本データvs一様分布（理想）のqqプロット n = ', num2str(n)])
end

%以下関数
function itiyo = itiyo(n) %ベルヌーイ分布から一様分布を作成する関数、nはループ回数
    sum = 0;
    
    for i = 1:n
        ber = binornd(1, 0.5);
        sum = sum + ber / (2^i);    
    end

    itiyo = sum;
end


function kai2_kentei = kai2_kentei(sample_size, ransu) %一様分布に従っているか、カイ二乗適合度検定のための統計量を算出する関数
    %初期値
    kai = 0;
    
    for i = 0:9
        expect = sample_size * 0.1; %期待度数
        obs = length(ransu( 0.1 * i < ransu & ransu < 0.1 * i + 0.1)); %観測度数
        kai = kai + (expect - obs)^2 / expect; %カイ二乗統計量を求める。
    end
    
    kai2_kentei = kai;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%box-muller法を使用し、標準正規乱数を生成するコード

n = 10

sample_size = 1000;

%一様乱数を生成する(samplesize ×2の配列）
for i = 1:2
    for j = 1:sample_size
        ransu(j,i) = itiyo(n);
    end
end

%box mullar法を用いて、一様乱数から正規乱数を生成する。
for s = 1 : sample_size
    X1(s) = sqrt( -2 * log(ransu(s,1))) * cos(2 * pi * ransu(s,2));
end

histogram(X1)
xlabel("乱数の値")
ylabel("度数")
title(['乱数の分布 n = ', num2str(n)], 'FontSize',12)

%qqプロットを行う
pd = makedist('Normal');

figure;
qqplot(X1,pd);

%コルモゴロフ-スミルノフ検定
h = kstest(X1)

if h == 0
    ans = "正規分布に従っていないとは言えない"
else
    ans = "正規分布に従っていない"
end

%以下関数
function itiyo = itiyo(n) %ベルヌーイ分布から一様分布を作成する関数、nはループ回数
    sum = 0;
    
    for i = 1:n
        ber = binornd(1, 0.5);
        sum = sum + ber / (2^i);    
    end

    itiyo = sum;
end

*1:可算無限回コインを投げられる人は各自で実験してみて下さい。

*2:もし乱数が一様分布に従っているのであれば、プロットされた点はほぼ45度線上に分布する。詳しくは以下のURLを参照 ja.wikipedia.org

*3:帰無仮説は「度数が（離散）一様分布に従っている」

*4:ぺこぱ

*5:あくまで離散化したデータに対してそう言えるだけで、元の連続的な乱数値に対しては何も言えないのに注意（階級の幅の取り方によって結果は変わる）

*6: (0,1)上の一様分布に従う確率変数X,Yを用いて、 \begin{align} Z_1 =\sqrt{-2 logX} cos 2 \pi Y \\ Z_2 =\sqrt{-2 logX} sin 2 \pi Y \\ \end{align} とすれば、Z_1とZ_2は独立に標準正規分布に従うというもの。 詳しくは ja.wikipedia.org

*7:ベクトルデータが標準正規分布に従っているか調べる検定。帰無仮説は「データが標準正規分布から発生している。」 jp.mathworks.com