Matlab蒙特卡罗模拟:
可以用蒙特卡罗方法来近似模拟求pi值:
思路:
设相互独立的随机变量X,Y均服从[-1,1]上的均匀分布,则(X,Y)服从{-1≤x≤1, 1≤y≤1}上的二元均匀分布(即图1中正方形区域上的二元均匀分布),记事件A = {x2+y2≤1},则A事件发生的概率等于单位圆面积除以边长为2的正方形的面积,即P(A) = pi/4,从而可得圆周率pi = 4P(A). 而P(A)可以通过蒙特卡洛模拟法求得,在图1中正方形内随机投点(即横坐标X和纵坐标Y都是[-1,1]上均匀分布的随机数),落在单位圆内的点的个数m与点的总数n的比值m/n可以作为A事件的概率P(A)的近似,随着投点总数的增加,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(‘投点个数’);
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; % 蒙特卡洛的次数
R = zeros(times,1); % 用来储存扰动项u和x1的相关系数
K = zeros(times,1); % 用来储存遗漏了x2之后,只用y对x1回归得到的回归系数
for i = 1: times
n = 30; % 样本数据量为n
x1 = -10+rand(n,1)*20; % x1在-10和10上均匀分布,大小为30*1
u1 = normrnd(0,5,n,1) - rand(n,1); % 随机生成一组随机数
x2 = 0.3*x1 + u1; % x2与x1的相关性不确定, 因为我们设定了x2要加上u1这个随机数
u = normrnd(0,1,n,1); % 扰动项u服从标准正态分布
y = 0.5 + 2 * x1 + 5 * x2 + u ; % 构造y
k = (n*sum(x1.*y)-sum(x1)*sum(y))/(n*sum(x1.*x1)-sum(x1)*sum(x1)); % y = k*x1+b 回归估计出来的k
K(i) = k;
u = x2 + u; % 因为我们忽略了x2,所以扰动项要加上x2
r = corrcoef(x1,u); % 2*2的相关系数矩阵
R(i) = r(2,1);
end
plot(R,K,'*')
xlabel("x_1和u'的相关系数")
ylabel("k的估计值")
结果: