読者です 読者をやめる 読者になる 読者になる

標準正規分布に従う乱数

perl 数値計算

標準正規分布に従う乱数を発生させるプログラムを書いてみました。
正規分布は次のような式になります。
f(x;\mu,\sigma) = \frac{1}{\sigma\sqrt{2\pi}}\exp{(-\frac{(x - \mu)^2}{2\sigma})}


標準正規分布は、「\mu=0, \sigma=1」の正規分布になります。
f(x; 0,1) = \frac{1}{sqrt{2\pi}}\exp{(-\frac{x^2}{2})}


この分布に従う乱数を生成するプログラムを書いて、生成結果を確かめてみました。

#!/usr/bin/env perl                                                                                
use strict;
use warnings;
use constant PI => 3.1415926535;


my $c = 1000000;


# 乱数生成と集計                                                                                   
my @r;
my %counter = ();
for my $i (0..$c) {
    my $f = gauss();
    $counter{int($f*10)/10}++;
    push @r, $f;
}


# 乱数の平均を求める                                                                               
my $sum = 0;
for my $r (@r) {
    $sum += $r;
}

# 乱数の分散を求める                                                                               
my $sum2 = 0;
for my $r (@r) {
    $sum2 += ($r - $ave) * ($r - $ave);
}
my $var = $sum2 / $c;
print "$var\n";


# ヒストグラムを表示する                                                                           
for my $m (sort { $a <=> $b } keys %counter) {
    my $count = $counter{$m} / 1000;
    my $kome = '*' x $count;
    if ($kome) {
        print $kome, "\n";
    }
}


# 標準正規分布を発生させる、逆関数
sub gauss {
    my $func = sqrt(-2.0*log(rand()))*cos(2*PI*rand());
    return $func;
}


出力結果は次のようになりました。

*
*
**
**
***
***
****
*****
*******
********
**********
***********
*************
***************
******************
********************
**********************
*************************
***************************
******************************
********************************
**********************************
***********************************
*************************************
***************************************
***************************************
*******************************************************************************
***************************************
**************************************
*************************************
************************************
**********************************
********************************
******************************
***************************
*************************
***********************
********************
******************
****************
**************
************
**********
********
*******
******
****
***
***
**
*
*
*


期待通りの分布になりました。