C言語による正規分布の乱数作成

C言語の勉強を始めたてで正規分布の乱数を作るプログラムを以下のプログラムを参考に勉強しているのですが、プログラムの内容や正規分布の公式がどの変数に当てはまるのか理解が追いついていません。ご教授願えれば幸いです。

C言語
1#include <stdio.h>
2#include <stdlib.h>
3#include <math.h>   
4
5double nrand();
6
7double nrand(){
8    
9    static int sw=0;
10    static double r1,r2,s;
11    
12    if (sw==0){
13        sw=1;
14        do {
15            r1=2.0*drand48()-1.0;
16            r2=2.0*drand48()-1.0;
17            s=r1*r1+r2*r2;
18        } while (s>1.0 || s==0.0);
19        s=sqrt(-2.0*log(s)/s);
20        return(r1*s);
21    }
22    else {
23        sw=0;
24        return(r2*s);
25    }
26}
27
28
29int main(){
30    
31    int i;
32    int seed = 1234567; 
33    double x;
34    
35    srand48( seed );
36    
37    for (i=0; i<1000; i++){ x = nrand();
38        printf("%.5lf\n",x);
39    }
40}
41

行動規範の内容に同意します

回答2件

これは、Box-Muller法という方法による正規分布乱数の生成です。

主役は完全に数学の世界の話で、プログラムは数式を起こしただけなので、詳細についてはWikipediaや数学サイトの記事などをご覧ください。

投稿2018/10/10 05:04

maisumakun

総合スコア145183

ベストアンサー

質問欄に記載のコードは Marsaglia polar 法です。
標準ライブラリ libstdc++ の std::normal_distribution もこのアルゴリズムで生成されています。

質問欄のコードの仕組み

-1 < x < 1, -1 < y < 1 の正方形からランダムに点 x, y を選びます。
これは、x, y が [-1, 1] の一様分布に従う観測値としたとき、s = x^2 + y^2 < 1 を満たす x, y を選ぶことを意味します。
これをやっているのが、コードのこの部分です。

c
1do {
2    r1 = 2. * drand48() - 1.;
3    r2 = 2. * drand48() - 1.;
4    s = r1 * r1 + r2 * r2;
5} while (s > 1. || s == 0.);

drand48() は [0, 1) の一様分布に従う乱数を返す関数のようです。
drand48() は C言語の標準APIではないので、rand() を使う場合、

c
1r1 = (double)rand() / RAND_MAX * 2. - 1.;
2r2 = (double)rand() / RAND_MAX * 2. - 1.;

となります。rand() は [0, RAND_MAX] の一様分布に従う乱数を返す関数なので、[0, RAND_MAX] -> [0, 1] -> [0, 2] -> [-1, 1] と範囲を変換しています。

このとき、s = r1^2 + r^2 とすると、標準正規分布に従う乱数を2つ次の式で生成できます。

証明はこちらにあります。

1回呼び出すと2つの乱数が生成できるので、1回呼び出したときは
r1 * sqrt(-2. * log(s) / s) を返して
次に呼び出したときは前回返してないもう一方の
r2 * sqrt(-2. * log(s) / s)
を返してます。
これらは static 変数なので、次に呼び出したときも値を参照できます。
この切り替えは sw をいう static 変数をフラグとして使っています。

if (sw == 0) {  // 生成したとき
    sw = 1;
    return r1 * s;
}
else {  // 次に呼び出したとき
    sw = 0;
    return r2 * s;
}

1000 個値を生成してプロットしてみた結果
仮説検定とかはやってませんが、正規分布になっていそうですね。

他のアルゴリズム

正規分布を生成するアルゴリズムは他にもあるようです。
いろいろなアルゴリズムの実装が提供されている Github を見つけたので、見てみると面白いかもしれません。

追記グラフの生成方法について

以下のコードでテキストファイルに出力して、Python で読み込んで seaborn というライブラリでプロットしました。
Cだけでプロットする場合は gunplot とか使うことになるのでしょうけど、それでも結構面倒くさいです。Python だと1行できれいなグラフが出力できます。

cpp
1#include <fstream>
2#include <math.h>
3#include <stdio.h>
4#include <stdlib.h>
5
6/**
7   @brief [-1, 1] の一様分布に従う乱数を生成する。
8   @return [-1, 1] の一様分布に従う乱数
9 */
10double urand()
11{
12    return (double)rand() / RAND_MAX * 2. - 1.;
13}
14
15double nrand()
16{
17    // 1度に2個の正規分布に従う値 r1, r2 が得られるので
18    // 生成したときは r1、次に呼び出したときは r2 を返す。
19    static int sw = 0;  // switch 用のフラグ
20    static double r1, r2, s;
21
22    if (sw == 0) {
23        sw = 1;
24        // (x, y) \in [-1, 1]^2 を満たす点をランダムに選ぶ。
25        do {
26            r1 = urand();
27            r2 = urand();
28            s = r1 * r1 + r2 * r2;
29        } while (s > 1.0 || s == 0.0);
30
31        // 正規分布に従う点に変換する。
32        s = sqrt(-2. * log(s) / s);
33        return r1 * s;
34    }
35    else {
36        sw = 0;
37        return r2 * s;
38    }
39}
40
41int main()
42{
43    int seed = 1234567;
44
45    srand(seed);
46    std::ofstream ofs("values.csv");
47    for (int i = 0; i < 1000; i++)
48        ofs << nrand() << std::endl;
49}

ファイルを読み込んで、出力する。

python
1import seaborn as sns
2import matplotlib.pyplot as plt
3import numpy as np
4
5array = np.loadtxt('values.csv')
6print(array.shape)
7
8sns.set()
9ax = sns.distplot(array)
10plt.show()