定積分 乱数

[ykneu]

[1] 質問文
　[1.1] 自分が今行いたい事は何か：モンテカルロ法で定積分∫(1,0) (1/sqrt(2π))*exp(-x^2/2)dx の近似値を計算したい。
　[1.4] 今何がわからないのか、知りたいのか：どのようにx、yを＜＝１にできるか分かりません。

[2] 環境 　
　[2.2] コンパイラ名 ： gcc

[box]

　[1.1] 自分が今行いたい事は何か：モンテカルロ法で定積分∫(1,0) (1/sqrt(2π))*exp(-x^2/2)dx の近似値を計算したい。
　[1.4] 今何がわからないのか、知りたいのか：どのようにx、yを＜＝１にできるか分かりません。

x = 0のときの被積分関数の値は1 / √(2π) ≒ 0.39894228
で、x = 1のときの被積分関数の値はe^(-1) / √(2π) ≒ 0.14676266
ですから、y <= 1というしばりはマズくないでしょうか。
yは、0から正規分布曲線上の範囲であればよい、ということにはならないでしょうか。

なお、[0, 1]の範囲の乱数を求める最も安直な方法は
標準関数のrand()を使って求めた0～RAND_MAXの値を
(double) RAND_MAX
で割ってやる、というものがあるように思います。
例えばメルセンヌ・ツイスターまで持ち出すのはちと大げさな気がします。

[box]

おもしろそうな問題だったので、ちょいとコードを書いてみました。

コード:

include Math

def f(x)
  exp(-x * x / 2) / sqrt(2 * PI)
end

N = 1000000
count = 0
N.times do
  x, y = rand, rand
  count += 1 if y <= f(x)
end
puts count.to_f / N

実行結果は、だいたい0.341前後です。
真の値は、手元の関数電卓によると、0.3413447461になりました。

[ykneu]

確かにyは正規分布表のうちであれば良いですね。
ご指摘ありがとうございました。
ただ、書いていただいたプログラムがよく理解できません。
くわしくおしえていただけませんか？

[かずま]

ただ、書いていただいたプログラムがよく理解できません。
くわしくおしえていただけませんか？

これなら分かりますか？

コード:

#include <stdio.h>   // printf
#include <stdlib.h>  // rand
#include <math.h>    // exp, sqrt

double inv_sqrt_2pi;

double f(double x) { return inv_sqrt_2pi * exp(x * x / -2); }

int main(void)
{
	int i, n = 1000000, count = 0;
	inv_sqrt_2pi = 1 / sqrt(2 * 3.141592653589793238);
	for (i = 0; i < n; i++) {
		double x = rand() / (RAND_MAX + 1.0);
		double y = rand() / (RAND_MAX + 1.0);
		if (y <= f(x)) count++;
	}
	printf("%f\n", (double)count / n);
	return 0;
}

x = 0のときの被積分関数の値は1 / √(2π) ≒ 0.39894228
で、x = 1のときの被積分関数の値はe^(-1) / √(2π) ≒ 0.14676266
ですから、y <= 1というしばりはマズくないでしょうか。

x = 1のときの被積分関数の値は、e^(-1/2) / √(2π) ≒ 0.24197072 ですよね。

y <= 1 というのは、0 <= x <= 1, 0 <= y <= 1 で、1 x 1 の正方形の面積だから、
count / n で値が求まるのではありませんか？

[box]

x = 1のときの被積分関数の値は、e^(-1/2) / √(2π) ≒ 0.24197072 ですよね。

おっと、間違えてしまいました。