Matlab蒙特卡羅模擬:
可以用蒙特卡羅方法來近似模擬求pi值:
思路:
設(shè)相互獨立的隨機(jī)變量X，Y均服從[-1,1]上的均勻分布，則(X,Y)服從{-1≤x≤1, 1≤y≤1}上的二元均勻分布（即圖1中正方形區(qū)域上的二元均勻分布），記事件A = {x2+y2≤1}，則A事件發(fā)生的概率等于單位圓面積除以邊長為2的正方形的面積，即P(A) = pi/4，從而可得圓周率pi = 4P(A). 而P(A)可以通過蒙特卡洛模擬法求得，在圖1中正方形內(nèi)隨機(jī)投點（即橫坐標(biāo)X和縱坐標(biāo)Y都是[-1,1]上均勻分布的隨機(jī)數(shù)），落在單位圓內(nèi)的點的個數(shù)m與點的總數(shù)n的比值m/n可以作為A事件的概率P(A)的近似，隨著投點總數(shù)的增加，m/n會越來越接近于P(A)，從而可以得到逐漸接近于pi的模擬值。
function piva=PiMonterCarlo(n)
x=0;y=0;d=0;
m=length(n);
pivalue=zeros(m,1);
for i=1:m
x=2rand(n(i),1)-1;
y=2rand(n(i),1)-1;
d=x.^2+y.2;
pivalue(i)=4sum(d<=1)/n(i);
end
if nargout==0
if m>1
figure;
plot(n,pivalue,‘k.’);
h=refline(0,pi);
set(h,‘linewidth’,2,‘color’,‘r’);
text(1.05n(end),pi,’\pi’,‘FontSize’,15);
xlabel(‘投點個數(shù)’);
ylabel(’\pi的模擬值’);
else
figure;
plot(x,y,‘r.’);
hold on;
h=rectangle(‘Position’,[-1 -1 2 2],‘LineWidth’,2);
t=linspace(0,2*pi,100);
plot(cos(t),sin(t),‘k’,‘LineWidth’,2);
xlabel(‘X’);ylabel(‘Y’);
title(‘Pi的模擬值:’,num2str(pivalue));
axis([-1.1 1.1 -1.1 1.1]);axis equal;
end
else
piva=pivalue;
end

還可以模擬K值

times = 300; % 蒙特卡洛的次數(shù) R = zeros(times,1); % 用來儲存擾動項u和x1的相關(guān)系數(shù) K = zeros(times,1); % 用來儲存遺漏了x2之后，只用y對x1回歸得到的回歸系數(shù) for i = 1: timesn = 30; % 樣本數(shù)據(jù)量為nx1 = -10+rand(n,1)*20; % x1在-10和10上均勻分布，大小為30*1u1 = normrnd(0,5,n,1) - rand(n,1); % 隨機(jī)生成一組隨機(jī)數(shù)x2 = 0.3*x1 + u1; % x2與x1的相關(guān)性不確定， 因為我們設(shè)定了x2要加上u1這個隨機(jī)數(shù)u = normrnd(0,1,n,1); % 擾動項u服從標(biāo)準(zhǔn)正態(tài)分布y = 0.5 + 2 * x1 + 5 * x2 + u ; % 構(gòu)造yk = (n*sum(x1.*y)-sum(x1)*sum(y))/(n*sum(x1.*x1)-sum(x1)*sum(x1)); % y = k*x1+b 回歸估計出來的kK(i) = k;u = x2 + u; % 因為我們忽略了x2，所以擾動項要加上x2r = corrcoef(x1,u); % 2*2的相關(guān)系數(shù)矩陣R(i) = r(2,1); end plot(R,K,'*') xlabel("x_1和u'的相關(guān)系數(shù)") ylabel("k的估計值")

結(jié)果:

總結(jié)

以上是生活随笔為你收集整理的Matlab蒙特卡罗模拟的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。