• ベストアンサー

モンテカルロ法

I=∫ f(x) dx (0から1までの定積分）を評価するために、一様分布U(0,1)に従うU1,...,Unを用いてIn=1/n(f(U1)+...+f(Un))を計算した。 a)E[In],V[In]を求めよ。 b)f(x)=xe^xとして、Inの相対誤差が5%以下になる確率を90%とするには、nをいくらにすればよいか。 ヒント(In-E[In]/sqrt(V[In]) ～ N(0,1)としてよい） とあって、非常に基本的な問題で申し訳ないのですが、a)がよくわかりません。意味的にn->∞の時、E[In]=Iであることはわかるのですが、a)はE[In]=Iとしてよいのでしょうか？

質問者が選んだベストアンサー

  • ベストアンサー
  • keyguy
  • ベストアンサー率28% (135/469)
回答No.5

V=∫(0<u<1)・(f(u))^2・duとおく Uの確率密度をp(u)とすると f(U)の平均は ∫(-∞<u<∞)・f(u)・p(u)・du =∫(0<u<1)・f(u)・du=I f(U)の分散は ∫(-∞<u<∞)・(f(u)-I)^2・p(u)・du =V-I^2 同様にして 確率変数U1,U2,U3,・・・,Unがそれそれ独立で密度がそれぞれp(u)とすると Inの平均はIとなり Inの分散は(V-I^2)/nとなる 確率変数(U1+U2+U3+・・・+Un)/nの密度は nが大きければ大きいほどN(I,(V-I^2)/n)に近い しかしこのような近似をしなくてもパソコンを使えば 厳密に密度関数を求める事ができる 正規分布で求める場合には この確率変数を正規化してN(0,1)を使う 0.1・I<|In-I|なる確率を0.05以下にすれば良い 正規化変数はZ=(In-I)/√((V-I^2)/n)であるから 0.1・I/√((V-I^2)/n)<|z|なる確率が0.05以下になるnを求めれば良い 正規分布表から両側の面積が0.05となる境を求めて その境=0.1・I/√((V-I^2)/n) とすればよい VとIを数値計算で求めておけば求まる

betagamma
質問者

お礼

遅くなって申し訳ありません。完璧な回答ありがとうございます！ 自分で計算したところ、全く同じ結果になりました。結局、n>=646ぐらいになるようです。

その他の回答 (5)

  • keyguy
  • ベストアンサー率28% (135/469)
回答No.6

VとIを数値計算で求めておけば求まる は VとIを求めておけばでる に変更 解析的に簡単にもとまる関数なので数値計算は不要ということ

betagamma
質問者

お礼

お礼が大変遅くなって申し訳ありません。 皆様、どうもありがとうございました。 そうですね、問題が数値計算をする問題なのに、目的の値まですべて解析的にもとまってしまう（練習用のためなのでしょうが）のでちょっと不安になりました。 レポートを提出した後、ほかの人のレポートも見せてもらったところ、やはり、自分が求めた値と同じぐらいになっていました。

  • yaksa
  • ベストアンサー率42% (84/197)
回答No.4

直接計算だと、ルベーグ積分が必要になる気がするけど、その気持ちだけ書くなら、 E[f(U)] = Σy P(y<f(U)<y+dy) = Σy P(U∈{x|y<f(x)<y+dy}) = Σy ρ({x|y<f(x)<y+dy}) = ∫f(x)dx = I ってことです。積分を縦の短冊の和として考えるんではなくて、横の短冊の和として考えるということ。 Ukは独立だから、E[In]=I 分散は平均の式で、f(x)の代わりに、g(x) = (f(x)-I)^2 を代入して、 V[f(U)] = E[g(U)] = ∫g(x)dx = ∫(f(x)-I)^2 dx なので、 V[In] = 1/n * ∫(f(x)-I)^2 dx ですね。

回答No.3

No1の回答で標準偏差となっているところは分散に訂正してください。

回答No.2

下の回答は書きかけで投稿してしまいました。「つまりいつもSnで」の後に書こうとしたのは次のようなことです。「Snで母集団の分散の推定とすると平均として真の分散より小さい推定値になります。」

回答No.1

Inは不偏な推定量なのでE[In]=Iとしてよいと思います。不偏でない推定量の代表は標準偏差 　Sn=Σ( Xn - Xbar )^2/n です。不偏標準偏差は 　S'n=Σ( Xn - Xbar )^2/(n-1) で定義されることは良く知られています。ｎを固定し、ｎ個の標本からSnを計算する操作をｋ回繰り返したとします。第i回目のSnをSniとして平均をSkとします。 　Sk=(Sn1+Sn2+…＋Snk)/k ｋ→∞としたとき、Skは母集団の分散には収束しません。つまりいつもSnで不偏でない推定量とはこのようなものだと思います。一方、Inについては、上のSkに相当するものを作り、(nでなく)ｋを無限大とするとIに収束します。私は統計の専門家ではなく、間違っているかもしれませんのでテキストブックで確認をお願いします。

  1. カテゴリ
  2. 学問・教育
  3. 数学・算数

関連するQ&A

  • ｆ（ｘ）＝1/sqrt(2*π*v))*e^(-(x-u)^2/(2*

    ｆ（ｘ）＝1/sqrt(2*π*v))*e^(-(x-u)^2/(2*v)について、ｆ’（ｘ）＝０、ｆ’’（ｘ）＝０となる点を求めよ。 という問題なのですが、答えはｆ’（ｘ）＝０はｘ＝u、ｆ’’（ｘ）＝０はｘ＝u-sqrt（ｖ）、u+sqrt（ｖ）でよいのでしょうか？ 答えを求めることを点を求めると解釈してよいのでしょうか？

  • シンプソン法の誤差評価

    (例)台形公式の誤差評価 区間[x(i-1) , x(i)]を考える 真の面積は、ΔS(i-1)=∫(x(i-1) →x(i)　）f(x)dx=F(x(i))-F(x(i-1)),　近似面積は、h/2(f(x(i-1)+f(x(i))である。 ところでF(x(i))=F(x(i-1)+h)=F(x(i-1))+f(x(i-1))h+(1/2!)f ' (x(i-1))h^2+,,,,,,であるから 誤差E(i-1)=真の面積-近似面積=ΔS(i-1)-f(x(i-1))h=(-1/12)f '' (x(i-1)h^3+,,,,, 全区間の誤差 E=Σ(i=1,n)E(i-1)=-n((h^3/12)Σ(i=1,n)f ' (x(i-1))+,,,,)=( (x(n)-x(0))/12 )h^2Σ(i=1,n)f ''(x(i-1)+,,,,,, と同様にしてシンプソンの公式の誤差評価も示したいのですが、解答、解法お願いします

  • 連立方程式をSOR法を用いてコンピュータで解く

    { 2 1 1 }{ X1 }　{ 2 } { 2 3 1 }{ X2 }＝{ 4 } { 1 1 3 }{ X3 } {-1 } という問題をSOR法を用いてコンピュータ（fortran）で解け。ただし初期値は０とし許容誤差eps=10^-5とする。そこで私は次のようなプログラムで解きました。 ____dimension a(20,20),b(20),x(20),x0(20) ____open(5,file='senkei.d') ____open(6,file='senkei.r') ____open(7,file='senkei.k') ____read(5,*)n,max,eps,u ____read(5,*)((a(i,j),i=1,n),j=1,n) ____read(5,*)(b(i),i=1,n) ____read(5,*)(x0(i),i=1,n) ____do i=1,n ____ad=a(i,i) ____a(i,i)=0.0 ____b(i)=b(i)/ad ____do j=1,n ____a(i,j)=a(i,j)/ad ____end do ____end do ____do k=1,max ____do i=1,n ____w=0.0 ____do j=1,n ____if(j.lt.i)x0(j)=x(j) ____w=w+a(i,j)*x0(j) ____end do ____x(i)=u*(b(i)-w)+(1-u)*x(i) ____write(7,10)x(i) _10_format(3f10.5) ____end do ____w=0.0 ____do i=1,n ____w=w+(x(i)-x0(i))**2 ____end do ____w=sqrt(w/n) ____if(w.lt.eps)go to 20 ____do i=1,n ____x0(i)=x(i) ____end do ____end do ____write(6,*)'error' ____go to 30 _20_write(6,*)'sulution' ____write(6,40)(x(i),i=1,n) _40_format(3f15.5) _30_stop ____end max:最大繰り返し回数,u:緩和係数,a(i,j):係数行列, b(i):既知項,x(i):未知項としました。データではeps=0.00001,max=50,緩和係数uは0.8から0.1ずつ増やして1.5までそれぞれ答えを求めました。私のプログラムではu=1.1の時7回で収束しました。しかしu=1.5の時3回で収束したのですが答えはかなりの誤差ありました。許容誤差eps=10^-5としたのになぜ誤差が大きいまま、しかも三回で収束してしまったのかわかりません。誰かわかる方いらっしゃったら教えていただけませんか？長々と申し訳ありません。お願いします。

  • 多次元正規分布に関する質問

    多次元正規分布についてのシミュレーションをする課題を持っているのですが 少しわからないところがあります。 m次元正規分布の式はwikipediaにあるように http://ja.wikipedia.org/wiki/%E6%AD%A3%E8%A6%8F%E5%88%86%E5%B8%83 N(u,S)=exp(-(x-u)'invS(x-u)/2) / (2PI)^(m/2)sqrt(detS) です。ですが分散共分散行列Sを使わず分散をスカラーとしてつまり σで書きたいのです。 （なのでm次元のどの方向においても分散はσ^2でいいです。 そういう簡単な状況でのシミュレーションをしようとしています。） この場合、m次元正規分布のしきはどうなるのでしょうか。 よろしくお願いいたします。

  • シンプソン法の出力結果について

    シンプソン法の区間分割数nを10～100まで10ずつ増やして計算値と計算誤差を求めるプログラムを書きに作成したのですが、出力結果に-7.341865999405491E-5などとあります。この「E-数字」とはJavaではどういうこと示しているのですか？ また、計算誤差の求め方は下記のプログラムでいいのですか？ public class Simpson { static double f(double x) { // ここに任意の被積分関数を記述 double y = Math.exp(- x * x / 2) / Math.sqrt(2.0 * Math.PI); return y; } public static void main(String[] args) { double a = - 1.0, b = 1.0; // 積分範囲 int n = 100; // 区間分割数 for(n=10; n <=100; n+=10){ double h = (b - a) / (double)n; // 分割幅 double s, s1 = 0.0, s2 = 0.0; for (int i = 1; i <= n / 2; i++) { s1 += f(a + (2 * i - 1) * h); } for (int i = 1; i <= n / 2 - 1; i++) { s2 += f(a + 2 * i * h); } s = h / 3.0 * (f(a) + 4.0 * s1 + 2.0 * s2 + f(b)); double suti = (s-0.68269)/0.68269*100;　　//計算誤差＝（計算値ー真値）/真　　 　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　値×100 System.out.println("区間分割数　＝" + n); System.out.println("シンプソン法による計算値　＝" + s); System.out.println("シンプソン法による計算誤差　＝"+suti+"\n"); } } }

  • 統計学でこまってます。

    レポートが不合格で返ってきました。どうしても わかりません。おしえてください。　　 　 離散型確率変数Ｘ、Ｙの分布は 　　Ｐ（Ｘ＝xi)=pi,P(Y=yi)=qi (i=1,2)です。 （１）Ｅ（Ｘ＋Ｙ）＝Ｅ（Ｘ）＋Ｅ（Ｙ） （２）ＸとＹが独立な確率変数であるとき 　　Ｖ（Ｘ＋Ｙ）＝Ｖ（Ｘ）＋Ｖ（Ｙ） 　批評は（１）Pi=ri1+ri2,qj=v1j+v2jを証明してください。 　　　　（２）Ｘ、Ｙが独立のとき 　　　　　　　Ｅ（ＸＹ）＝Ｅ（Ｘ）Ｅ（Ｙ）を証明する。 確率変数Ｘが二項分布Ｂ（９、１/2)に従う時、Xの分布の値 Ｐ（Ｘ＝k)=(0～９）のひとつひとつを正規分布で近似し 相対誤差を計算する。ここで相対誤差｜ｄ/P(X=k)|*100%, d:誤差です。数値は小数点以下第6位を四捨五入して第5位まで。 批評はＰ（Ｘ＝Ｋ）[K=0～９]のひとつひとつを正規近似する。 　　　9Ｃ0=(1/2)^0=1です。 上記3問よろしくお願いします。

  • モンテカルロ法の面積近似

    Mが2の32乗の、線形合同法で乱数を発生し、 モンテカルロ法を使って、 領域a = {(x,y)|x≧1,y≧2,(x-1)^2 + (y-2)^2 ≦1}の面積の近似値を計算するプログラムを作成し、これを利用してn=10000とした場合の上の領域の面積、および円周率πの近似値を求めよ。 (関数を使う) という課題が出たのですが、乱数の出し方？が変であまり近似されません(^_^;) どなたかどこが変か教えてください(>_<) これだと、2と4になってしまうんです(/_;) #include<stdio.h> #include<stdlib.h> #include<math.h> #define N 100 /*領域Rを含む面積T　領域R=求める面積*/ #define Nx 1 #define Ny 2 /*乱数を発生させる関数*/ unsigned long int random_g(unsigned long int x0, unsigned long int a, unsigned long int c); int main(){ unsigned long int x0, a, c, M; double x, y, s; int n, count_in, count_out; count_in = 0; count_out = 0; printf("x0: ", x0 ); scanf("%d", &x0); a = 61; c = 49; for(n = 0; n < 10000; n++){ /*乱数を発生*/ x = (random_g(x0, 4*a*n+1, c) / (double)4294967296)*Nx + 1; y = (random_g(x0, 4*a*n+1, c) / (double)4294967296)*Ny + 1; if(x >= 1 && y >= 1 && pow(x-1, 2) + pow(y-2, 2) <= 1) count_in++; else count_out++; } s = (double)count_in/(count_in + count_out)*Nx*Ny; printf("s = %f\n", s); printf("pi = %f\n", 2*s); return 0; } unsigned long int random_g(unsigned long int x0, unsigned long int a, unsigned long int c){ int i; x0 = (a*x0 + c); /*printf("%u\n", x0);*/ return x0; }

  • ２つの媒介変数がある関数が表す領域の図示

    uは区間[0:a](a>0)を、vは区間[0:1]を動くとします。ここで、f(u,v)=uv/sqrt(1+v^2) + v , g(u,v)=-u/sqrt(1+v^2) + v^2/2　として、領域Ω＝｛(x,y)|x=f(u,v), y=g(u,v)｝を図示せよという問題なんですが、どのように考えるのが定石でしょうか？（領域の形は分かっています）

  • 密度関数

    X，Yの同時分布は密度関数2exp(-x-y)，0＜x＜y＜∞(その他では0)をもつ。 (a)　このとき、X，Yのそれぞれの周辺分布の密度関数ｆ1(ｘ)，ｆ2(ｙ)および条件X=4のもとでのYの密度関数f(y|X=4)を求めよ。 (b)　同時分布U=X+Y，V=X/Yの密度関数g(u，v）を求めよ。U，Vは独立になるか。 (a)の方は自分で解けました。 (b)の方が分からないのでどなたか教えてください。 答えはg(u，v)=2{e^(-u)}u/(1+v)^2 でU,Vは独立になるようです。

  • モンテカルロ法を用いた実験

    モンテカルロ法を用いてy=x^2とx軸で囲まれた部分の面積の近似値を求める．これを利用して平均μ，標準偏差σ，を求めたところ下記の写真のようになった．点数nに対して標準偏差は何次であるかを，中心極限定理から求められる次数，実験結果から導かれる次数のそれぞれで求めなさい．次数とはどのように求められるのでしょうか？使用したプログラムを下に書いておきます． #include <stdio.h> #include <time.h> #include <stdlib.h> #include <math.h> double f(double x); int main(void) { int M = 1000; int i,j,n,c,k; int tmp[6] = { 20, 50, 100, 200, 500, 1000 }; double x,y,Y, sumY, sumY2, muY, muY2, sigmaY; srand((unsigned int)time(NULL)); for (k = 0; k < 6; k++){ sumY = 0; sumY2 = 0; for (j = 0; j < M; j++){ c = 0; n = tmp[k]; for (i = 0; i < n; i++){ x = (1.0* (double)rand()) / ((double)RAND_MAX + 1.0); y = (1.0* (double)rand()) / ((double)RAND_MAX + 1.0); if (y <= f(x)){ c++; } } Y = ((double)c) / (double)n; sumY = sumY + Y; sumY2 += Y*Y; } muY = sumY / (double)M; muY2 =sumY2 / (double)M; sigmaY = sqrt(muY2 -(muY*muY)); printf("%d %lf %lf\n", n, muY, sigmaY); } return 0; } double f(double x) { double r = x * x; return r; }

注目のQ&A

カテゴリ

専門家に質問してみよう

職業から探して質問する

あなたにピッタリな商品が見つかる！ OKWAVEセレクト

質問する