正規乱数の生成

中心極限定理

中心極限定理とは下記のとおり。

赤本からの引用です。

統計学入門 (基礎統計学Ⅰ)

統計学入門 (基礎統計学Ⅰ)

  • 発売日: 1991/07/09
  • メディア: 単行本


母集団分布が何であっても、確率変数X_1,X_2,\cdots,X_nの和S=X_1+X_2+\cdots+X_nの確率分布の形はnが大なるときには大略正規分布と考えてよい。

例えば、上記の確率変数Xが平均\mu、分散\sigma^2に従うとする時、

S=X_1+X_2+\cdots+X_n \sim N(n\mu,n\sigma^2)

もしくは、

\bar{X}=S/n \sim N(\mu,\sigma^2/n)


Sの期待値と分散は次のようになります。\bar{X}についても同様に導出できます。

E[S] = E[X_1+X_2+\cdots+X_n] =E[X_1]+E[X_2]+\cdots+E[X_n] =n\mu

V[S] = V[X_1]+V[X_2]+\cdots+V[X_n] = n\sigma^2

なんとなく、どんな母集団の分布でも標本数を増やすと正規分布になる、みたいな誤解がありそうですが、そんなことはなく、誤解を恐れずに言うと、標本の和や平均が正規分布に近づいていくということです。

普通に考えて、母集団が対数正規分布に従うのに、その母集団からサンプリングすると正規分布になるなんてありえません。

標準正規分布に従う乱数の生成

ボックスミューラーの正規乱数生成については紹介しましたが、今回は上記の中心極限定理を使って生成したいと思います。

中心極限定理では、元の母集団はなんでもよいとのことですので、一様分布に従う確率変数X\sim U(0,1)を考えます。この一様分布の確率変数Xn個サンプリングすると、その平均\bar{X}N(\mu,\sigma^2/n)に従います。

次にXは一様分布U(0,1)に従うので、その期待値と分散は次のようになります。

nishiru3.hatenablog.com

E[X] = (a+b)/2 = (0+1)/2 = 1/2
V[X] = (b-a)^2/12 = (1-0)^2/12 = 1/12

以上から、一様分布からn個のサンプルを拾ってきて和を取ったS正規分布N(n/2,n/12)に従います。

さらにこれを標準化します。実際の値をrとすると、その和は\sum_{i=1}^{n} r_i と表現できるので、標準化した値zは次のようになります。

z=\dfrac{\sum_{i=1}^{n} r_i - n/2}{\sqrt{n/12}}

標準化したので、zは標準正規分布N(0,1)に従います。では実装してみましょう。

#!/usr/bin/env perl
use strict;
use warnings;

# perl random.pl > random_out.dat

my $n = 12;
my @random_numbers;

for my $i ( 1..1000 ) { # 1000個の乱数を生成
    my $sum_r = 0.0; 
    for my $j ( 1..$n ) {
        $sum_r += rand(1); # 一様乱数の生成と足し合わせ
    }
    my $z = ($sum_r - $n/2)/sqrt($n/12);
    push @random_numbers, $z;
}
for my $k( @random_numbers ) {
    print $k."\n";
}

gnuplotを使ってヒストグラムを書きました。引用も追加しておきます。

set terminal png
set output 'random_out.png'
binwidth =0.5
bin(x,width)=width*floor(x/width)+width/2.0
plot [-4:4] 'random_out.dat' using (bin($1,binwidth)):(1.0) smooth freq with boxes title 'random\_out.dat'

nucl.phys.s.u-tokyo.ac.jp

f:id:nishiru3:20200229225701p:plain

それっぽく生成されています。
ちなみに、今回はnを12として、1000個の正規乱数を生成しています。実はn=12とすると、z=\dfrac{\sum_{i=1}^{n} r_i - n/2}{\sqrt{n/12}}=\sum_{i=1}^{n} r_i - 6になるので、式が簡単になり、平方根の計算は入りません。コードは簡略化していませんが。

変数変換による正規分布N(\mu,\sigma^2)に従う乱数の生成

標準正規分布N(0,1)の乱数を作りましたが、次にN(\mu,\sigma^2)のような一般的な正規分布を作ります。

と言っても非常に簡単です。変数変換を施して生成します。生成したい平均\mu、分散\sigma^2の正規乱数をx、標準正規分布に従う変数をzとします。
z=\dfrac{x-\mu}{\sigma}
x = \mu + \sigma z

上記のzを生成して、\sigmaを乗じて、\muを足し込めば完成です。さて実装です。

#!/usr/bin/env perl
use strict;
use warnings;

# perl random_g.pl > random_out_g.dat

my $n = 12;
my $mu = 10.0;
my $sigma2 = 5.0;
my $sigma = sqrt($sigma2);

my @random_numbers;

for my $i ( 1..1000 ) { # 1000個の乱数を生成
    my $sum_r = 0.0;
    for my $j ( 1..$n ) {
        $sum_r += rand(1); # 一様乱数の生成と足し合わせ
    }
    my $z = ($sum_r - $n/2)/sqrt($n/12);
    my $x = $mu + $sigma * $z;
    push @random_numbers, $x;
}
for my $k( @random_numbers ) {
    print $k."\n";
}
set terminal png
set output 'random_out.png'
binwidth =0.5
bin(x,width)=width*floor(x/width)+width/2.0
plot [0:20] 'random_out_g.dat' using (bin($1,binwidth)):(1.0) smooth freq with boxes title 'random\_out\_g.dat'

f:id:nishiru3:20200229232653p:plain

最後に

ここに書いているアルゴリズムを使って実際に利用するのはやめた方が良いです。ライブラリを使いましょう。PerlならPDLに関数ran_gaussianというものがあるようです。