rand()関数による一様乱数

    モンテカルロ法のシミュレーション等で一様乱数やその他の分布を持つ乱数が必要となることがあります。ここではまず、[0,1] 間の一様乱数を生成する方法について説明します。
    一番簡単なのが、C言語の組み込み関数を使う方法でしょう。rand() 関数を用いると、0 以上 RAND_MAX 以下の疑似乱数を発生させることができます。よって、rand()+0.5 を RAND_MAX+1 で割れば [0,1] の一様乱数になります（rand() を RAND_MAX で割るのではなく、0.5 を足して RAND_MAX+1 で割ることにより、両端に偏らず均等に分布するように出来ます）。
    しかし、RAND_MAX の値は処理系によっては比較的小さい（例えば32767）ので、そのままでは取り得る値の間隔が粗すぎて良くありません。そこで、以下のように、複数の rand() 関数の値を用いて一様乱数を計算すれば間隔の細かい乱数を生成することができます。
    #include <stdio.h>
#include <stdlib.h>
#include <time.h>

double urand(){
    double m, a;
    m = RAND_MAX + 1;
    a = (rand() + 0.5)/m;
    a = (rand() + a)/m;
    return (rand() + a)/m;
}

int main(){
    int i;
    srand((unsigned int)time(NULL));
    for(i=0; i<10; i++){
        printf("%f\n", urand());
    }
    return 0;
}
    RAND_MAX=32767=215-1 の場合、上のプログラムで呼び出している rand() 関数の3つの値が独立なら、2-45 の分解能の乱数が実現できます。実際は、rand() 関数の周期が通常 231 しかないため、245通りの数値は発生しません。
    ただ、rand()関数による疑似乱数は、統計的にはあまり良い性質を持っていない、環境により実装が異なり再現性に難がある、等の問題があるため、モンテカルロ法のシミュレーションには適しているとは言えません。

メルセンヌツイスタによる一様乱数

    メルセンヌツイスタは、松本眞氏と西村拓士氏によって開発された、高品質の疑似乱数を生成するアルゴリズムです。モンテカルロ法に用いるのに十分な品質を持っており、219937-1 という長い周期を持つことが証明されています。http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/mt.html に解説やソースコード等があります。
    メルセンヌツイスタのアルゴリズムはいくつかありますが、整数を返すアルゴリズムの他にも、倍精度浮動小数点数を返す高速なアルゴリズムもありますので、そちらを用いることで直接、一様乱数を生成することができます。また、速度を気にしなければ、整数版の結果を 上の rand() 関数を用いる場合のように実数化しても良いでしょう。

Xorshiftによる一様乱数

    Xorshiftは、2003年に George Marsaglia によって提案された疑似乱数発生アルゴリズムで、整数型変数のビット操作だけで性質の良い疑似乱数を高速に計算することができます。Xorshiftも、モンテカルロ法に用いるには十分な性能を持っています。ただ、メルセンヌツイスタと比較すると品質は若干見劣りするようです。
    Xorshiftには複数のバージョンがありますが、以下は周期が 2128-1 のものです。メルセンヌツイスタに比べると周期は短いですが、それでも、1秒に1京回の演算を行うスーパーコンピュータで計算しても周期が一巡するのに7500兆年ほどかかりますので、普通に使う分には十分に長い周期であると言えます。0 以上 232-1 以下の整数が返ります。
    #include <stdint.h>

uint32_t xorshift(){
    static uint32_t x = 123456789;
    static uint32_t y = 362436069;
    static uint32_t z = 521288629;
    static uint32_t w = 88675123;
    uint32_t t;
    t = x ^ (x<<11);
    x = y; y = z; z = w;
    w ^= t ^ (t>>8) ^ (w>>19);
    return w;
}
    これを実数化し、[0,1] の一様乱数を返すようにしたものが以下になります。上のプログラムでは乱数の系列は固定でしたが、下のプログラムでは seed を与えて系列を選べるようにしています。
    #include <stdint.h>
#include <stdio.h>
#include <time.h>

static uint32_t xorsft_x = 123456789;
static uint32_t xorsft_y = 362436069;
static uint32_t xorsft_z = 521288629;
static uint32_t xorsft_w = 88675123;

void initrand(uint32_t seed){
    do {
        seed = seed*1812433253 + 1; seed ^= seed<<13; seed ^= seed>>17;
        xorsft_x = 123464980 ^ seed;
        seed = seed*1812433253 + 1; seed ^= seed<<13; seed ^= seed>>17;
        xorsft_y = 3447902351 ^ seed;
        seed = seed*1812433253 + 1; seed ^= seed<<13; seed ^= seed>>17;
        xorsft_z = 2859490775 ^ seed;
        seed = seed*1812433253 + 1; seed ^= seed<<13; seed ^= seed>>17;
        xorsft_w = 47621719 ^ seed;
    } while(xorsft_x==0 && xorsft_y==0 && xorsft_z==0 && xorsft_w==0);
}

double urand(){
    uint32_t t;
    t = xorsft_x ^ (xorsft_x<<11);
    xorsft_x = xorsft_y;
    xorsft_y = xorsft_z;
    xorsft_z = xorsft_w;
    xorsft_w ^= t ^ (t>>8) ^ (xorsft_w>>19);
    return ((xorsft_x+0.5) / 4294967296.0 + xorsft_w) / 4294967296.0;
}

int main(){
    int i;
    initrand((unsigned int)time(NULL));
    for(i=0; i<10; i++){
        printf("%f\n", urand());
    }
    return 0;
}
    乱数の系列を初期化する initrand() では、近い値の seed を与えても出来るだけ独立な系列が得られるように工夫してあります。また、初期状態は initrand(0) と同じ状態になっています。

正規乱数の発生

    平均 0、分散 1 の正規分布 N[0,1] に従う正規乱数を発生させたい場合、基本的には [0,1] の一様乱数から生成することになります。よく用いられるのが Box-Muller 法です。X と Y を [0,1] の一様分布に従う確率変数とすると、Z1=√(-2log X) * cos(2πY) と Z2=√(-2log X) * sin(2πY) が N[0,1] に従う互いに独立な確率変数となります。
    プログラムとしては以下のようになります。Box-Muller 法では、互いに独立な正規乱数のペアを求めることができますが、ここでは片方のみ用いています。また、urand() 周りのプログラムは省略していますので、上で示したコードを適宜追加して下さい。
    #include <math.h>
#include <stdio.h>

double normrand(){
    return sqrt(-2*log(urand()))*cos(2*M_PI*urand());
}

int main(){
    int i;
    initrand((unsigned int)time(NULL));
    for(i=0; i<10; i++){
        printf("%f\n", normrand());
    }
    return 0;
}
    以下は、Marsaglia polar 法と言われる方法で、三角関数を計算しない分、Box-Muller 法よりも高速です（の筈ですが、私が試したところ Box-Muller 法と速度は大差ありませんでした）。
    #include <math.h>
#include <stdio.h>

double normrand(){
    double x, y, r;
    do{
        x = 2*urand() - 1;
        y = 2*urand() - 1;
        r = x*x + y*y;
    } while(r>=1);
    return sqrt(-2*log(r)/r)*x;
}

int main(){
    int i;
    initrand((unsigned int)time(NULL));
    for(i=0; i<10; i++){
        printf("%f\n", normrand());
    }
    return 0;
}

指数分布に従う乱数の発生

    パラメータ λ の指数分布に従う乱数を発生させたい場合は、逆関数法を用います。具体的には、X を [0,1] の一様分布に従う確率変数とすると、Z=-λlog(X) が パラメータ λ の指数分布に従う確率変数となります。
    プログラムは以下のようになります。正規乱数のときと同様、urand() 周りのプログラムは省略していますので、適宜、コードを追加して下さい。
    #include <math.h>
#include <stdio.h>

double exprand(double lambda){
	return -lambda*log(urand());
}

int main(){
    int i;
    initrand((unsigned int)time(NULL));
    for(i=0; i<10; i++){
        printf("%f\n", exprand(1.0));
    }
    return 0;
}

ポアソン分布に従う乱数の発生

    パラメータ λ のポアソン分布に従う乱数を発生させたい場合は、[0,1] の一様分布に従う k 個の互いに独立な確率変数の積が exp(-λ) を下回る確率が λkexp(-λ)/k! となることを利用します。
    プログラムは以下のようになります。正規乱数のときと同様、urand() 周りのプログラムは省略していますので、適宜、コードを追加して下さい。
    #include <math.h>
#include <stdio.h>

int porand(double lambda){
    double a, b;
    int k = 0;
    a = urand();
    b = exp(-lambda);
    while(a>=b){
        a *= urand();
        k++;
    }
    return k;
}

int main(){
    int i;
    initrand((unsigned int)time(NULL));
    for(i=0; i<10; i++){
        printf("%d\n", porand(1.0));
    }
    return 0;
}

全部合わせたプログラム

    以上で説明したコードを全部まとめたプログラムを以下に示します。ただし、一様乱数は Xorshift をベースとしています。
    #include <stdint.h>
#include <stdio.h>
#include <math.h>
#include <time.h>

static uint32_t xorsft_x = 123456789;
static uint32_t xorsft_y = 362436069;
static uint32_t xorsft_z = 521288629;
static uint32_t xorsft_w = 88675123;

// 乱数系列初期化
void initrand(uint32_t seed){
    do {
        seed = seed*1812433253 + 1; seed ^= seed<<13; seed ^= seed>>17;
        xorsft_x = 123464980 ^ seed;
        seed = seed*1812433253 + 1; seed ^= seed<<13; seed ^= seed>>17;
        xorsft_y = 3447902351 ^ seed;
        seed = seed*1812433253 + 1; seed ^= seed<<13; seed ^= seed>>17;
        xorsft_z = 2859490775 ^ seed;
        seed = seed*1812433253 + 1; seed ^= seed<<13; seed ^= seed>>17;
        xorsft_w = 47621719 ^ seed;
    } while(xorsft_x==0 && xorsft_y==0 && xorsft_z==0 && xorsft_w==0);
}

// 一様乱数
double urand(){
    uint32_t t;
    t = xorsft_x ^ (xorsft_x<<11);
    xorsft_x = xorsft_y;
    xorsft_y = xorsft_z;
    xorsft_z = xorsft_w;
    xorsft_w ^= t ^ (t>>8) ^ (xorsft_w>>19);
    return ((xorsft_x+0.5) / 4294967296.0 + xorsft_w) / 4294967296.0;
}

// 正規乱数
double normrand(){
    return sqrt(-2*log(urand()))*cos(2*M_PI*urand());
}

// 指数分布に従う乱数
double exprand(double lambda){
	return -lambda*log(urand());
}

// ポアソン分布に従う乱数
int porand(double lambda){
    double a, b;
    int k = 0;
    a = urand();
    b = exp(-lambda);
    while(a>=b){
        a *= urand();
        k++;
    }
    return k;
}

int main(){
    int i;
    initrand((unsigned int)time(NULL));
    printf("一様乱数\n");
    for(i=0; i<10; i++){
        printf("%f\n", urand());
    }
    printf("正規乱数\n");
    for(i=0; i<10; i++){
        printf("%f\n", normrand());
    }
    printf("指数分布に従う乱数\n");
    for(i=0; i<10; i++){
        printf("%f\n", exprand(1.0));
    }
    printf("ポアソン分布に従う乱数\n");
    for(i=0; i<10; i++){
        printf("%d\n", porand(1.0));
    }
    return 0;
}