定積分 乱数

フォーラム(掲示板)ルール
フォーラム(掲示板)ルールはこちら  ※コードを貼り付ける場合は [code][/code] で囲って下さい。詳しくはこちら
ykneu

定積分 乱数

#1

投稿記事 by ykneu » 8年前

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

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

box
記事: 2002
登録日時: 13年前

Re: 定積分 乱数

#2

投稿記事 by box » 8年前

ykneu さんが書きました:  [1.1] 自分が今行いたい事は何か:モンテカルロ法で定積分∫(1,0) (1/sqrt(2π))*exp(-x^2/2)dx の近似値を計算したい。
 [1.4] 今何がわからないのか、知りたいのか:どのようにx、yを<=1にできるか分かりません。
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
記事: 2002
登録日時: 13年前

Re: 定積分 乱数

#3

投稿記事 by box » 8年前

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

コード:

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

Re: 定積分 乱数

#4

投稿記事 by ykneu » 8年前

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

かずま

Re: 定積分 乱数

#5

投稿記事 by かずま » 8年前

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 で値が求まるのではありませんか?

box
記事: 2002
登録日時: 13年前

Re: 定積分 乱数

#6

投稿記事 by box » 8年前

かずま さんが書きました: x = 1のときの被積分関数の値は、e^(-1/2) / √(2π) ≒ 0.24197072 ですよね。
おっと、間違えてしまいました。
バグのないプログラムはない。
プログラムは思ったとおりには動かない。書いたとおりに動く。

閉鎖

“C言語何でも質問掲示板” へ戻る