累積平均の求め方を知りたいです

ビュフォンの針という、円周率をシュミレーションによって求める問題を研究しているmatlab初心者の学生です。

得られた値の累積平均が実際の円周率の値にどのように収束するか調べたいのですが、どのように累積平均を取ればいいのかわかりません。

教えていただきたいです。

N=1000;%針の数
L=0.20;%針の長さ
xb=L+rand(1,N)*(1-2*L);%針の始点のx座標
yb=L*rand(1,N)*(1-2*L);%針の始点のy座標
angs=rand(1,N)*360;%針の傾き
xe=xb+L*cosd(angs);
ye=yb+L*sind(angs);
ax=axes;
plot(ax,[xb;xe],[yb;ye]);%針を表示
axis square
hold on
glines=0:L:1;%平行線を表示
for i =1:length(glines)
    xline(ax,glines(i));
end

n=sum(floor(xb/L)~=floor(xe/L));%平行線に交わった針の数

n = 627

piEstimate=2*N/n%円周率の値を求める

piEstimate = 3.1898

毎回変数を変え、このスクリプトで得られたpiEstimateの値を、実行するたびに今までの値累積平均をとって保存したいと思っています。

方法を教えていただけるとありがたいです。

Answer 1

[Atsushi Ueno]

piEstAry = []; % 最初に空ベクトルとして初期化しておく
for i = 1:5
    piEstAry = [piEstAry, i] % 過去の演算値をベクトルに追記していく
    piEstimate = mean(piEstAry) % 平均値を得る
end

piEstAry = 1

piEstimate = 1

piEstAry = 1×2

1 2

piEstimate = 1.5000

piEstAry = 1×3

1 2 3

piEstimate = 2

piEstAry = 1×4

1 2 3 4

piEstimate = 2.5000

piEstAry = 1×5

1 2 3 4 5

piEstimate = 3

[健]

ありがとうございました!

無事得られた結果の累積平均を求めることができました。

[Atsushi Ueno]

過去に演算した値を残す必要は無いと思います。その場合は下記の方法が良いと思います。なお、for文のイテレータを i にすると複素数と混同するので k に変更しました。

piEstimate = 0; % 最初に0で初期化しておく
...
for k = 1:5 % 繰り返し回数 = 累積数
    piEstimate = ((k-1) * piEstimate + k) / k % 平均値を得る
end

piEstimate = 1

piEstimate = 1.5000

piEstimate = 2

piEstimate = 2.5000

piEstimate = 3

[健]

ありがとうございます！私の理解不足だと思うのですが、教えていただいたコードで実行したところ、 結果が表示されませんでした。averagepiのところが怪しいと思っているのですが、どうすればよいでしょうか？ コードは以下になります。よろしくお願いします

if true
  % code
k=10;%試行回数を設定
for s=1:k
N = 1000*s;%針の数
L = 0.20;%針の長さ
xb = L + rand(1,N)*(1-2*L);%針の始点のx座標
yb = L + rand(1,N)*(1-2*L);%針の始点のy座標
angs = rand(1,N)*360;%針の傾き
xe = xb + L*cosd(angs);%判定用の値
ye = yb + L*sind(angs);
ax = axes;
plot(ax,[xb;xe],[yb;ye]);%針を表示
axis square
hold on
glines = 0:L:1;%平行線を表示
for i = 1:length(glines)
xline(ax, glines(i));
end
n = sum(floor(xb/L) ~= floor(xe/L));%平行線に交わった針を求める
piEstimate = 2 * N / n;%円周率の値を求める
averagepi=0;
averagepi=((k-1)*piEstimate+k)/k;

end

[Atsushi Ueno]

回答の「過去に求めた円周率リスト」と、後からコメントした「過去に演算した値を残さずに求める累積平均」を両方とも計算しています。プロットすると時間が掛かるのでコメントアウトしています。

試行回数を1000回にして「過去に求めた円周率リスト」をプロットすると下記の様になりました。最後に求めたaveragepi = 3.141511994945709 と期待通りになりました。

L = 0.20;%針の長さ
k = 1000;%試行回数を設定
averagepi = 0; % 最初に0で初期化しておく
piEstAry = []; % 最初に空ベクトルとして初期化しておく
for s = 1:k
    N = 1000*s;%針の数
    xb = L + rand(1,N)*(1-2*L);%針の始点のx座標
    yb = L + rand(1,N)*(1-2*L);%針の始点のy座標
    angs = rand(1,N)*360;%針の傾き
    xe = xb + L*cosd(angs);%判定用の値
    ye = yb + L*sind(angs);
%     ax = axes;
%     plot(ax,[xb;xe],[yb;ye]);%針を表示
%     axis square
%     hold on
%     glines = 0:L:1;%平行線を表示
%     for i = 1:length(glines)
%         xline(ax, glines(i));
%     end
    n = sum(floor(xb/L) ~= floor(xe/L));%平行線に交わった針を求める
    piEstimate = 2 * N / n;%円周率の値を求める
    piEstAry = [piEstAry, piEstimate]; % 過去の演算値をベクトルに追記していく
    averagepi = ((s - 1) * averagepi + piEstimate) / s;%求めた円周率の累積平均を求める
end

[Atsushi Ueno]

averagepiを求める為にpiEstAryは不要である事に注意ください。

[健]

そうですね！ありがとうございました。