MENU

正規乱数表現

正規乱数

ja.wikipedia.org

m を平均, \sigma を標準偏差とすると, 正規分布 N(m,\sigma^{2}) は次のように表現される分布となります。

{ \displaystyle N(m,\sigma^{2}) = \frac{1}{\sqrt{2 \pi \sigma^{2}}} exp^{-\frac{(x-m)^2}{2 \sigma^{2}}} }

f:id:Manao55:20201228230325p:plain
画像引用元：「Metrology: The Science of Measurement」https://www.muelaner.com/wp-content/uploads/2013/07/Standard_deviation_diagram.png

この平均 m , 標準偏差 \sigma の正規分布 N(m,\sigma^{2}) に従う乱数は, ボックス・ミュラー法により, 2個の一様乱数 r_1 , r_2 を用いて

{ \displaystyle x_1 = \sigma \sqrt{-2 logr_1} cos(2 \pi r_2) + m }
{ \displaystyle x_1 = \sigma \sqrt{-2 logr_1} sin(2 \pi r_2) + m }

で得ることができます。

プログラム

#include <iostream>  /* cout */
#include <math.h>    /* sin, cos */
#include <cstdlib>   /* rand */
#include <iomanip>   /* setw */
using namespace std;

void boxrnd(double, double, double *, double *);

int main() {
    int i, j, hist[1000];
    double x,y;

    for (i=0;i<1000;i++)
        hist[i] = 0;

    for (i=0;i<10000;i++) {
        boxrnd(2.6,10.0,&x,&y);
        hist[(int)x]++;
        hist[(int)y]++;
    }

    for (i=0;i<20;i++){
        cout << setw(3) << i+1 << " : " << setw(4) << hist[i] << " I " ;
        for(j=1;j<=hist[i]/100;j++){
            cout << "*" ;
        }
        cout << endl;
    }
    return 0;
}

void boxrnd(double sig, double m, double *x, double *y){

    double r1,r2;
    r1 = rand()/32767.1;
    r2 = rand()/32767.1;
    *x = sig*sqrt(-2*log(r1))*cos(2*3.14159*r2) + m;
    *y = sig*sqrt(-2*log(r1))*sin(2*3.14159*r2) + m;
}

実行結果

  1 :    5 I
  2 :   17 I
  3 :   50 I
  4 :  141 I *
  5 :  335 I ***
  6 :  730 I *******
  7 : 1267 I ************
  8 : 1939 I *******************
  9 : 2613 I **************************
 10 : 3069 I ******************************
 11 : 2987 I *****************************
 12 : 2528 I *************************
 13 : 1897 I ******************
 14 : 1207 I ************
 15 :  669 I ******
 16 :  333 I ***
 17 :  140 I *
 18 :   47 I
 19 :   19 I
 20 :    6 I