コンピュータ上で「ランダムさ」を再現するには疑似乱数を用いる。本課題では、標準一様乱数 \(U \sim U(0,1)\) を起点として、用途に合う 4 つの確率分布に従う乱数を生成する。
下のテキストエリアには 配布/sample1.cpp の中身がそのまま貼り付けてある。// *** の箇所を埋めて、
「sample1.cpp をダウンロード」で自分の環境に保存できる。
配布/sample1.cpp）の // *** を埋める:
exp_dist は 逆関数法: \(y = -\dfrac{1}{\lambda}\ln(1-u)\)normal_dist は Box–Muller 法: \(\theta = 2\pi u_2,\; z = \sqrt{-2\ln(1-u_1)}\cos\theta\)lognormal_dist は パラメタ変換: \(v = \ln(1+\mathrm{var}_0/m_0^2),\; m = \ln m_0 - v/2\)DIST_TYPE を 1〜4 に変えて 4 回コンパイル・実行:
g++ sample1.cpp && ./a.out > data.txt末尾の
Average: ..., Variance: ... 行はツールが自動スキップするのでそのまま使える。
data.txt を読み込み、度数分布表と理論 PDF を比較する。メインループでは 1 行 1 サンプルを標準出力に、末尾で平均・分散のまとめ行を出力する。
0.289067 0.195331 0.741802 ... Average: 0.500151, Variance: 0.083120
Average: 行は HTML ツールが自動で除外するので、そのまま渡せる。
g++ sample1.cpp ./a.out > data.txt
data.txt を指定、分布とパラメタ（sample1.cpp での設定と同じ値）を入れて「集計」を押す。
度数分布表・ヒストグラム棒・理論 PDF 曲線が重なって表示される。
DIST_TYPE | 1:一様, 2:指数, 3:正規, 4:対数正規。4 つとも動作させる必要あり。 |
NUM_TRIAL | 生成するサンプル数（既定 100,000）。大きいほど度数分布が滑らかに。 |
seed | 線形合同法の初期値 \(x_0\)。再現性のため固定（例: 19）。 |
PI | Box–Muller で使用する \(\pi\)。 |
漸化式 \( x_{n+1} = (a \cdot x_n + c) \bmod M \) で整数列を作り、\(M{-}1\) で割って \([0,1)\) の値に変換する。
double uniform_dist(long *x){
long M=65536, a=997, c=1;
*x = (a*(*x)+c) % M;
return (double)*x/(M-1);
}
学術的な品質は高くない（周期 \(\le 65536\)）が、仕組みの学習用。長時間シミュレーションでは sample2/3 の ran2()（周期 \(\approx 2\times 10^{18}\)）を使う。
指数分布の累積分布関数 \( F(y) = 1 - e^{-\lambda y} \) を \(y\) について解いて \( y = F^{-1}(u) \) を得る:
\[ y = -\frac{1}{\lambda} \ln(1 - u), \quad u \in [0, 1) \]u==1.0 の場合は \(\log(0)\) が発散するので do-while で排除している。
double exp_dist(long *x, double lam){
double u;
do { u = uniform_dist(x); } while (u == 1.0);
double y = 0;
// *** y の右辺を指数乱数が発生するように変形
return y;
}
2 つの独立な一様乱数 \(u_1, u_2 \in [0,1)\) から標準正規乱数 \(z \sim \mathcal{N}(0,1)\) を次で得る:
\[ \theta = 2\pi\, u_2, \qquad z = \sqrt{-2 \ln(1-u_1)} \cos\theta \]平均 \(m\), 分散 \(\mathrm{var}\) の正規乱数は \( \sqrt{\mathrm{var}}\, z + m \)。return sqrt(var)*z+m の部分は埋まっている。
double normal_dist(long *x, double m, double var){
double z, u1, u2, theta;
do { u1 = uniform_dist(x); } while (u1 == 1.0);
do { u2 = uniform_dist(x); } while (u2 == 1.0);
theta = 0;
// *** theta の右辺を変形
z = 0;
// *** z の右辺を変形
return sqrt(var)*z + m;
}
平均 \(m_0\), 分散 \(\mathrm{var}_0\) の対数正規分布を得るには、内部の正規分布のパラメタを次で設定する:
\[ \mathrm{var} = \ln\!\left(1 + \tfrac{\mathrm{var}_0}{m_0^2}\right), \qquad m = \ln m_0 - \tfrac{\mathrm{var}}{2} \]その正規乱数 \(z\) を exp(z) すれば、目的の対数正規乱数が得られる。
double lognormal_dist(long *x, double m0, double var0){
double m = 0, var = 0;
//*** 対数正規分布の平均 m0, 分散 var0 となるように
//*** 正規分布の平均 m, 分散 var を設定
double z = normal_dist(x, m, var);
return exp(z);
}
| 分布 | パラメタ | 確率密度関数 \(f(x)\) | 平均 | 分散 |
|---|---|---|---|---|
| 一様 | \(\alpha, \beta\) | \( \dfrac{1}{\beta-\alpha} \quad (\alpha \le x \lt \beta) \) | \( \tfrac{\alpha+\beta}{2} \) | \( \tfrac{(\beta-\alpha)^2}{12} \) |
| 指数 | \(\lambda\) | \( \lambda e^{-\lambda x} \) | \( 1/\lambda \) | \( 1/\lambda^2 \) |
| 正規 | \(m, \mathrm{var}\) | \( \tfrac{1}{\sqrt{2\pi\mathrm{var}}} \exp\!\big({-}\tfrac{(x-m)^2}{2\,\mathrm{var}}\big) \) | \( m \) | \( \mathrm{var} \) |
| 対数正規 | \(m_0, \mathrm{var}_0\) | \( \tfrac{1}{x\sqrt{2\pi v}} \exp\!\big({-}\tfrac{(\ln x - \mu)^2}{2v}\big) \) \(v = \ln(1{+}\mathrm{var}_0/m_0^2),\ \mu = \ln m_0 - v/2\) |
\( m_0 \) | \( \mathrm{var}_0 \) |
case 1: 標準一様乱数 \(u\) を \([\alpha, \beta]\) の一様乱数に変換して rnum に代入exp_dist: \( y = -\frac{1}{\lambda} \ln(1-u) \)normal_dist: \(\theta = 2\pi u_2\)normal_dist: \( z = \sqrt{-2\ln(1-u_1)} \cos\theta \)lognormal_dist: var, m を対数正規の与えられた \(m_0, \mathrm{var}_0\) から計算