ページ 1 ／ 1
定積分 乱数
Posted: 2015年11月12日(木) 20:26
by ykneu
[1] 質問文
[1.1] 自分が今行いたい事は何か：モンテカルロ法で定積分∫(1,0) (1/sqrt(2π))*exp(-x^2/2)dx の近似値を計算したい。
[1.4] 今何がわからないのか、知りたいのか：どのようにx、yを＜＝１にできるか分かりません。
[2] 環境
[2.2] コンパイラ名 ： gcc
Re: 定積分 乱数
Posted: 2015年11月12日(木) 22:54
by box
ykneu さんが書きました:
[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
で割ってやる、というものがあるように思います。
例えばメルセンヌ・ツイスターまで持ち出すのはちと大げさな気がします。
Re: 定積分 乱数
Posted: 2015年11月12日(木) 23:32
by 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になりました。
Re: 定積分 乱数
Posted: 2015年11月14日(土) 00:27
by ykneu
確かにyは正規分布表のうちであれば良いですね。
ご指摘ありがとうございました。
ただ、書いていただいたプログラムがよく理解できません。
くわしくおしえていただけませんか？
Re: 定積分 乱数
Posted: 2015年11月14日(土) 01:31
by かずま
ykneu さんが書きました:
ただ、書いていただいたプログラムがよく理解できません。
くわしくおしえていただけませんか？
これなら分かりますか？
コード:
#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;
}
box さんが書きました:
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 で値が求まるのではありませんか？
Re: 定積分 乱数
Posted: 2015年11月14日(土) 02:06
by box
かずま さんが書きました:
x = 1のときの被積分関数の値は、e^(-1/2) / √(2π) ≒ 0.24197072 ですよね。
おっと、間違えてしまいました。