Box-Muller法で発生させた乱数で正規分布を描く

[Under Construction]

Box-Muller法で正規分布に従う乱数を発生させて、その乱数で実際に正規分布をヒストグラムで描画するプログラムを作る。

これは、平均0、分散1の標準正規分布。

プログラムは、「distribution.c」

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

#define NT     10000000             // 発生させる乱数の数
#define BOX    500                  // ヒストグラムの解像度
#define UNIT   0.02                 // 刻み幅
#define WIDTH  10                   // WIDTH := BOX*UNIT
#define LEDGE -0.5*WIDTH            // 最小の座標(Left Edge)
#define VALUE  (double)BOX/NT/WIDTH // ヒストグラムの単位量

void   setSeed() {srand((unsigned)time(NULL));}
double genRand() {return rand()/(1.0+RAND_MAX);}

double gaussian(double myu, double sgm)
{
  static int    frag = 0;
  static double save = 0.0;
  
  if(frag == 0){
    double u_0 = genRand();
    double u_1 = genRand();
    double v_0 = myu + sgm*sqrt(-2*log(1-u_0))*cos(2*M_PI*u_1);
    double v_1 = myu + sgm*sqrt(-2*log(1-u_0))*sin(2*M_PI*u_1);
    save       = v_1;
    frag       = 1;
    return v_0;
  }
  else{
    frag = 0;
    return save;
  }  
}

void initArray(double array[], int size)
{
  int it;
  for(it = 0; it < size; it++)
    array[it] = 0.0;
}

void histogram(double hist[], int size, double nRand)
{
  double ibox      = (nRand-LEDGE)/UNIT;
  hist[(int)ibox] += VALUE;
}

double calcMean(const double nRand[], int size)
{
  double sum = 0.0;
  int    it;
  for(it = 0; it < size; it++)
    sum += nRand[it];
  return sum/size;
}

double calcVariance(const double nRand[], int size, double mean)
{
  double sum = 0.0;
  int    it;
  for(it = 0; it < size; it++)
    sum += nRand[it]*nRand[it] - mean*mean;
  return sum/size;
}

void plotter(FILE* gp, const double hist[], int size)
{
  int it;
  fprintf(gp, "set terminal x11\n");
  fprintf(gp, "set title 'Normal Distribution'\n");
  fprintf(gp, "unset key\n");
  fprintf(gp, "plot '-' w l\n");
  for(it = 0; it < size; it++){
    fprintf(gp, "%g %g\n", LEDGE+UNIT*it, hist[it]);
  }
  fprintf(gp, "end\n");
  fflush(gp);
}

int main(int argc, char *argv[])
{
  static double hist[BOX]; // ヒストグラム
  static double nRand[NT]; // 乱数
  int    it;
  double myu    = 0.0; // Box-Muller法で発生させる分布の平均
  double sgm    = 1.0; // Box-Muller法で発生させる分布の標準偏差
  double mean   = 0.0; // 実際に計算する分布の平均
  double vari   = 0.0; // 実際に計算する分布の分散
  FILE   *gp    = popen("gnuplot -persist", "w");
  
  if(2 <= argc && argc <= 3){
    myu = atof(argv[1]);
    sgm = atof(argv[2]);
    sgm = sqrt(sgm);
  }
  else if(argc >= 4){
    fprintf(stderr, "Error! Arguments too many!!\n");
    exit(EXIT_FAILURE);
  }
  setSeed();
  initArray(hist, BOX);
  initArray(nRand, NT);

  for(it = 0; it < NT; it++){
    nRand[it] = gaussian(myu, sgm);
    histogram(hist, BOX, nRand[it]);
  }

  mean = calcMean(nRand, NT);
  vari = calcVariance(nRand, NT, mean);

  plotter(gp, hist, BOX);
  printf("========================================\n");
  printf("Mean = %g, Variance = %g\n", mean, vari);
  printf("========================================\n");
  pclose(gp);

  return 0;
}

引数に平均と分散の値を入力して実行すると、それに従った分布が描画される。(Wikipedia：正規分布)

ヒストグラムを作る関数は工夫しないと実行速度が極端に遅くなる。