正規乱数の生成 - nishiru3の日記

中心極限定理

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

赤本からの引用です。

発売日: 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)$ を考えます。この一様分布の確率変数 $X$ を $n$ 個サンプリングすると、その平均 $\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'