乱数を使う数値計算のプログラム例

πの計算

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

#include "zdsfmt.c"  // 乱数関数の組み込み (プロジェクトに登録してあれば，必要ない）

int main() {

  int N = 100, i,j, NinCircle;
  int Nsample[5] = {100, 1000, 10000, 100000, 10000000};
  double x, y;

  InitMt(43431); // 乱数の種

  for(i=0; i<5; i++) {
    N = Nsample[i];
    n_in_circle = 0;  // 円の中に入った数

    for(j=0; j<N; j++) {
     x = NextUnif();
     y = NextUnif();
     if(x*x + y*y <= 1.0) n_in_circle++;
    }
    x = (double) n_in_circle/N * 4.0;
    printf("N=%d, pi=%lf, diff=%lf \n", N, x, x - acos(-1.));
  }
}

ランダムウォーク

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

#include "zdsfmt.c"  // 乱数関数の組み込み
#define WALK_COUNT 100 // １人が歩く歩数

void outDistribution(int *freq, int walkCount, int nPersons);
double average(int *freq, int walkCount, int nPersons);
double stdDiv(int *freq, int walkCount, int nPersons);
double normalDist(double x, double average, double variance);

int main() {
  int freq[WALK_COUNT*2+1] = {0}; // 度数記憶用配列
  int nPersons = 100000; // 試行する人数
  int position = 0;    // 歩行者の位置
  int i,j;

  InitMt(434231);

  for(i=0; i<nPersons; i++) {
    position = 0;
    for(j = 0; j<WALK_COUNT; j++) {
      if(NextUnif()<0.5) position++;
      else position--;
    }
    freq[position+WALK_COUNT]++;
  }

  //
  printf("平均=%lf, 標準偏差=%lf, サンプル数N=%d, sqrt(歩数)=%lf\n",
         average(freq, WALK_COUNT, nPersons), stdDiv(freq, WALK_COUNT, nPersons),
         nPersons, sqrt(WALK_COUNT));
  // 分布の出力
  outDistribution(freq, WALK_COUNT, nPersons);

}

/**** 度数分布の出力 *****************/
void outDistribution(int *freq, int walkCount, int nPersons) {
  int i;
  double m = 0.0, sigma = (double) walkCount; // 理論平均値と分散

  for(i=0; i<=walkCount*2; i+=2) {
    printf("%d, %lf, %lf\n", i-walkCount, (double)freq[i]/nPersons,
           2.0*normalDist((double)i-walkCount,m, sigma));
  }
}

/**** 平均値の計算 *******************/
double average(int *freq, int walkCount, int nPersons) {
  int i, a=0;
  for(i=0; i<=walkCount*2; i++) {
    a += freq[i]*(i-walkCount);
  }
  return (double) a/nPersons;
}

/**** 標準偏差の計算 *****************/
double stdDiv(int *freq, int walkCount, int nPersons) {

  int i, a2=0, a1=0;
  for(i=0; i<=walkCount*2; i++) {
    a1 += freq[i]*(i-walkCount);
    a2 += freq[i]*(i-walkCount)*(i-walkCount);
  }
  return sqrt((double) a2/nPersons - ((double)a1/nPersons));
}

/**** 正規分布 f(x), 引数 x, 平均値average，分散variance*****/
double normalDist(double x, double average, double variance) {
  double pi = acos(-1.0);
  return 1.0/(sqrt(2.0*pi*variance) * exp(x*x/2.0/variance));
}

乱数列の質の検証

上にあげた解説のページ

http://www.sat.t.u-tokyo.ac.jp/~omi/random_variables_generation.html

の2.1.4にしたがってrand()の検証をしてみましょう。

以下のように、rand()を [0,1] の乱数を発生させ、連続する３つの値を(x,y,z)として鳥瞰図に描いてみましょう。図示することで非一様性がよくわかります。 [0,1] のうち [0,0.01] の範囲に入った場合のみ出力します。（全域を大雑把にみれば一様に見えますが、この範囲では一様性が失われています。）

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

#define N 200
#define OUT_RANGE  0.01  // 出力する値の上限値

int main() {
  FILE *FP;
  int i=0;
  double x[3];

  FP = fopen("rand3d.csv","w");

  x[0]=(double)rand()/RAND_MAX;
  x[1]=(double)rand()/RAND_MAX;
  x[2]=(double)rand()/RAND_MAX;

  while(i < 1000) {
        if(x[0]<OUT_RANGE && x[1]<OUT_RANGE && x[2]<OUT_RANGE) {
          fprintf(FP, "%f, %f, %f\n",x[0],x[1],x[2]);
          i++;
        }
        x[0] = x[1];
        x[1] = x[2];
        x[2] = (double)rand()/RAND_MAX;
  }
  fclose(FP);
}

３次元の散布図を書く

scilabの場合

param3dを用います。

// csvで区切られた３次元座標データの点プロット
val = csvRead('z:\research\prog3ee\rand3d.csv') // csvファイル名(各自の設定)
param3d(val(:,1), val(:,2), val(:,3))

obj = gce()
obj.line_mode ="off"
obj.mark_style=0
obj.mark_size=1

matlabの場合

scatter3関数を用います。

val = csvread('z:\research\prog3ee\rand3d.csv')
scatter3(val(:,1),val(:,2),val(:,3),6,[1 0 0],'fill')

マウスドラッグにより視点を変えてみると非一様性がわかるはずです。