モンテカルロ法の積分教えてください。

[オス豚]

2重積分なんですけど、

∬f(x)dxdy
ただしf(x)=√(9-x^2-9/4*y^2)
　　　R={(x,y); 4x^2+9y^2≦36}

という問題です。Rというのは∬の右下に小さく書いてあるものです。
f(x)はなんとか自分でも作れるんですが、2重積分の計算ができません。

[管理人]

2重積分は、２重ループで計算できますよ。
普通の積分はとけますか？
小さくきって、四角形を作っていき、その合計を出す方法。
2重積分は、１方向をとめて、１方向ずつ、加算していきます。

ちょっと出かける用事がありますので、帰宅後、お答えします。
すみません。

[オス豚]

ある枠内に一様に乱数の点を発生させて、その図形の内部に入る割合によって面積を求める方法なんです。
普通の積分は何とかできます。

忙しい中すいません。

[管理人]

う、、モンテカルロ法昔やったんですけど、
すっかり忘れてしまいました・・。
いまや部分積分すら危うい予感・・；

調べてみて

http://momonga.t.u-tokyo.ac.jp/~ooura/intde-j.html

この辺が近いかな？と思うんですけど、どうでしょうか。
力不足でごめんなさい(_ _||)

[オス豚]

何とかできるかがんばって見ます。ありがとうございました。

[kuz]

初めまして、kuzといいます。
最近Cから離れてしまっていて、あっている自信は
ありません(というか間違えてる可能性が高いです)
しかもスマートじゃありません・・・・
ただ、お答えがない様ですので投稿いたします。
参考程度に見ていただければ幸いです

まず、Rの両辺を36で割ると
X^2/3^2 + Y^2/2^2 <=1 となるので
これは楕円の方程式をあらわします。
ここから-3<=X<=3,-2<=Y<=2
であることがわかります。
この範囲内の乱数の組み合わせで
Rを満たすものの組をf(x)に代入した結果の
総計を組数で割った値が近似値となります。

乱数を出力する所が見苦しいですが、以下ソースです

#include<stdio.h>
#include<stdlib.h>
#include<time.h>       
#include<math.h>       

#define TRY_CNT (10000) //最大試行回数

void main(void){

	double randX,randY,ans=0.0f;
	int i;


	for(i=0;i<TRY_CNT;i++){
		srand(time(NULL));   //乱数初期化
		for(;;){             //条件に合致するまで
			randX = (double)(rand() % 3001);
			randY = (double)(rand() % 2001);
			if(rand() % 2){
				randX = randX / 10000;
			}else{
				randX = randX / 100000;
			}
			if(rand() % 2){
				randY = randY / 10000;
			}else{
				randY = randY / 100000;
			}
			if(rand() % 2){
				randX = randX * (-1);
			}
			if(rand() % 2){
				randY = randY * (-1);
			}
			
			if(4*randX*randX + 9*randY*randY <= 36){
				break;
			}
		}
		
		ans += sqrt(9-randX*randX - 9/4*randY*randY);  //試行回数分集計
	}
	
	ans /= TRY_CNT;  //試行回数で割る
	
	printf("ans=%f",ans);
}

[box]

> randX = (double)(rand() % 3001);
> randY = (double)(rand() % 2001);

randXは0～3000、randYは0～2000の範囲の値をとります。

次の、10000または100000で割っている理由がよくわかりません。
いずれにせよ、randXは高々0.3、randYは高々0.2です。

その次に符号を変えるロジックがありますが、
後で2乗しているので不要だと思います。

> if(4*randX*randX + 9*randY*randY <= 36){
> break;

randXの2乗は高々0.09、randYの2乗は高々0.04です。
したがって、ここは必ずbreakします。
forループを回る回数は、常に1回です。

> ans += sqrt(9-randX*randX - 9/4*randY*randY); //試行回数分集計

ansの値は、3にひじょうに近くなります。
求める二重積分の値が楕円の面積だとすると、
π＊長径＊短径＝6π になるはずです。

[kuz]

すいません説明不足でした。

R(範囲)が楕円の方程式になっているというこです。
ここでRというのはx,yが満たすべき条件です。
そこで、本来なら適当に乱数を発生させ、
それがＲを満たすかどうかを判断し、
満たした場合f(x)に代入し、その総計
(TRY_CNT分)を取得します。そして集計を
TRY_CNTで除算した値が結果となります。

ただ、Rが楕円であることから発生させるべき
乱数の範囲がそれぞれ-3<x<3,-2<y<2ということなのです。
そこでその範囲内の乱数を発生させるために、
randXの場合では2000までの値を1000または10000で
除算することで0.001～2.000での乱数を発生させています。
Cで少数の乱数を発生させる方法が他に考えつかなかったので、
こんな、見苦しいことをしています、申し訳ありません

負についは確かに不要ですね。範囲内の乱数という意味で
負もあるほうが良いかと思ったのですが、蛇足でした。

まだ、わかりにくい点多々ありそうですが、
あくまで参考程度ということで、ご容赦ください

[オス豚]

ありがとうございます。
レベルの高い話しであまりついていけませんけども笑、参考にさせてもらってできました。