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的估计值")

结果:

	

蒙特卡罗算法，到底是干什么的？这个名词经常在金融行业或者其他许多领域听到，今天就和大家揭开这个神秘的面纱。
为什么会在很多领域听到蒙特卡罗，因为这是一个很神奇的算法，它可以模拟出很多场景，并且模拟出来的数据，可能与真实的数据相差无几，虽然不是十分精准，但是可以用，给人的感觉就是：
“我也不知道为什么，反正这样搞就能解决”
但这样模拟的成本，远远低于真实数据的获取的成本。或者说，模拟的这些数据，在实际生活中，暂时也是无法获取到的，比如对未来经济的预期。
废话不多少，今天，就用蒙特卡罗算法，做两个简单的模拟。一个是π值计算，另外一个求积分。
一、π值
π值是一个无理数，无限不循环，早在南北朝时期，我国数学家祖之冲得出精确到小数点后7位的结果。今天，我们用计算机来模拟一把，看看结果如何。
模拟思路：
如下图所示，可以推到正方形和内切圆的面积，存在比例关系，只要计算出它俩的面积比值，我们就可以求出π。我们可以用打点的方式，在正方形区域随机打点n个，如果落在内切圆的区域内有x个，则它俩的面积比就是n/x。如果这个n无限大时，则结果无限趋近于π。
 具体实现的核心代码如下：
while j < first_count: #first_count，模拟的总次数  x = random.uniform(-1, 1)  y = random.uniform(0, 1)  if   x**2 + y**2 :               count_s = count_s + 1 #在圆内的点的个数    j = j + 1
　　模拟出来的结果如下，在模拟超过1w次后，结果已经趋于稳定，基本等于3.14，这已经基本满足我们大部分使用场景了。
二、积分
积分实际也可以理解是计算区域内面积，比如下图，是y = -x^2+1 的函数图形，现在用蒙特卡罗求一下该函数的积分。
思路和求π的方法一致，也是通过随机打点的方式，根据在积分区域的散点数与矩形区域内散点数之比，乘以矩形面积，就是该积分区域面积。
分析模拟结果如下图，可以看到模拟3w到多次时，准确率很高了，与1.33不断接近，在9w次之后，基本保持重叠。
 通过蒙特卡罗模拟，生成一系列符合预期要求的随机数，就可以模拟出一个十分接近实际值的近似值，一般适应于对数值计算精度要求不是很高的场景，比如，我们在计算圆面积时，通常都会取3.14，而不会取3.1415926.....等。
☞ 推荐阅读 ☜
1、文本数据可以这样导入数据库
2、凯利公式助你获得更多收益
3、苹果最近的股价有点......
4、crontab定时任务了解一下
5、此av非彼"AV"
6、Oracle大规模数据快速导出文本文件
	

问题补充：
能提供一个例子看看吗？我就不懂MATLAB，想知道具体代码。因为现在写论文急，也没时间仔细看书了
蒙特卡罗模拟
就是随机数相关的东西，你只要知道随机数是怎么得到。其它的事就要好办了。
rand(m,n)产生m*n均匀随机数。
ex:
用概率方法求pi
N=100000;
x=rand(N,1);
y=rand(N,1);
count=0;
for i=1:N
if (x(i)^2+y(i)^2<=1)
count=count+1;
end
end
PI=4*count/N
-------------------------------
有关使用matlab进行蒙特卡罗模拟的程序问题
|
2010-10-1 15:37
|
浏览次数：1425次
已经分布是均匀分布(连续)，区间为(12,62)，请问各位大侠，如何用matlab编程实现此蒙特卡罗模拟，我想模拟2000次，得到概率密度图与累积概率密度图，程序应该如何编，麻烦大家指教，期待帮助，谢谢。小妹现在没有分数，还是希望大家能帮我，谢谢1
问题补充：
n=2000; %随机点数(可增加点数)
x=12+(62-12)*rand(1,n); %产生2000个12到62的随机数
xx=12:2:62; %画概率密度图的区间
nx=histc(x,xx);  %计算x在xx每个小区间内的点数。
px=nx/n;
sumpx=cumsum(px);
subplot(1,2,1)
bar(xx(1:end-1),px(1:end-1));
title('概率密度')
subplot(1,2,2)
plot(xx(1:end-1),sumpx(1:end-1));
title('累积概率密度')
-----------------------------------
我这儿正好有份程序，希望有所帮助。
代码：
一份蒙特卡洛程序
count=input('input the count:');%输入模拟粒子数
sigmaedata=[28370,13845,6908,2555,1223,2602,1925,905.3,479.7,164.5,74.24,23.86,66.60,36.62,22.29,9.978,5.298,1.668,0.7378,0.2361,0.1099,0.06211,0.03939,0.02030];
sigmacdata=[0.0220,0.0393,0.0568,0.0904,0.1209,0.1479,0.1722,0.2136,0.2480,0.3092,0.3486,0.3932,0.4153,0.4268,0.4319,0.4291,0.4215,0.3969,0.3691,0.3269,0.2944,0.2709,0.2512,0.2209];
Edata=[1,1.5,2,3,4,5,6,8,10,15,20,30,40,50,60,80,100,150,200,300,400,500,600,800];%截面数据
channel=zeros(1,ceil(662/5)+10);%多道数组
nget=0;%探测到的总计数
ntotal=0;%进入探测器的总计数
for ii=1:count%count个粒子循环
collidetime=0;%当前粒子碰撞次数
%粒子状态初始化
E0=622;
E=E0;
z=-2;
r=0;
theta=2*pi*rand(1);% 源抽样，z,r,theta坐标
miu=2*rand(1)-1;
fai=2*pi*rand(1);%方向角抽样
if miu
continue;
else
z=0;
r=2*sqrt(1-miu^2)/miu;
theta=fai;
end
while E>1%一个粒子在闪烁体中的输运过程
sigmae=interp1(Edata,sigmaedata,E,'linear');
sigmac=interp1(Edata,sigmacdata,E,'linear');
sigmat=sigmae+sigmac;%线性插值得到截面数据
L=-log(rand(1))/sigmat;%下次碰撞的距离
%计算下次碰撞位置坐标
rnew=sqrt(r^2+L^2*(1-miu^2)+2*r*L*sqrt(1-miu^2)*cos(fai-theta));
z=z+L*miu;
cdth=(rnew^2+r^2-L^2*(1-miu^2))/2/r/rnew;
sdth=L*sqrt(1-miu^2)*sin(fai-theta)/rnew;
dtheta=asin(sdth);
if cdth<0
dtheta=pi-dtheta;
end
theta=theta+dtheta;
r=rnew;
if(r>2)|(z>=4)|(z<0)%判断是否在闪烁体内
break;
else
collidetime=collidetime+1;
end
if rand(1)
E=0;
else%康普顿散射
alpha=E/511;
flag=0;
while flag==0
if rand(1)<=27/(4*alpha+29)
x=(1+2*alpha)/(1+2*alpha*rand(1));
if rand(1)<=0.5*((alpha+1-x/alpha)^2+1)
flag=1;
end
else
x=1+2*alpha*rand(1);
if rand(1)<=27/4*((x-1)^2)/x^3
flag=1;
end end
end E=E/x;
alphat=alpha/x;
miuL=1-1/alphat+1/alpha;%散射后的方向
a=miuL;
b=sqrt(1-a^2);
randangle=2*pi*rand(1);
miunew=a*miu+b*sqrt(1-miu^2)*cos(randangle);
sdf=b*sin(randangle)/sqrt(1-miunew^2);
cdf=(a-miu*miunew)/sqrt(1-miu^2)/sqrt(1-miunew^2);
sfn=sdf*cos(fai)+cdf*sin(fai);
cfn=cdf*cos(fai)-sdf*sin(fai);
fainew=asin(sfn);
if(cfn<0)

	

微信公众号：数据皮皮侠如果你觉得该公众号对你有帮助，欢迎关注、推广和宣传
内容目录：蒙特卡洛模拟方法及Python实现
1.什么是蒙特卡洛方法(Monte Carlo method)2.蒙特卡洛方法的基本思想3.蒙特卡洛求定积分4.蒙特卡洛方法python实例
1.什么是蒙特卡洛方法(Monte Carlo method)
       蒙特卡罗方法也称统计模拟方法，是1940年代中期由于科学技术的发展和电子计算机的发明，而提出的一种以概率统计理论为指导的数值计算方法。是指使用随机数(或更常见的伪随机数)来解决很多计算问题的方法。
      20世纪40年代，在冯·诺伊曼，斯塔尼斯拉夫·乌拉姆和尼古拉斯·梅特罗波利斯在洛斯阿拉莫斯国家实验室为核武器计划工作时，发明了蒙特卡罗方法。因为乌拉姆的叔叔经常在摩纳哥的蒙特卡洛赌场输钱得名，而蒙特卡罗方法正是以概率为基础的方法。
      与它对应的是确定性算法。
2.蒙特卡洛方法的基本思想
      通常蒙特卡罗方法可以粗略地分成两类：一类是所求解的问题本身具有内在的随机性，借助计算机的运算能力可以直接模拟这种随机的过程。例如在核物理研究中，分析中子在反应堆中的传输过程。中子与原子核作用受到量子力学规律的制约，人们只能知道它们相互作用发生的概率，却无法准确获得中子与原子核作用时的位置以及裂变产生的新中子的行进速率和方向。科学家依据其概率进行随机抽样得到裂变位置、速度和方向，这样模拟大量中子的行为后，经过统计就能获得中子传输的范围，作为反应堆设计的依据。
      另一种类型是所求解问题可以转化为某种随机分布的特征数，比如随机事件出现的概率，或者随机变量的期望值。通过随机抽样的方法，以随机事件出现的频率估计其概率，或者以抽样的数字特征估算随机变量的数字特征，并将其作为问题的解。这种方法多用于求解复杂的多维积分问题。
3.蒙特卡洛求定积分
蒙特卡洛方法的一个重要应用就是求定积分。来看下面的一个例子(参考文献2)
     当我们在[a,b]之间随机取一点x时，它对应的函数值就是f(x)。接下来我们就可以用f(x) * (b - a)来粗略估计曲线下方的面积，也就是我们需要求的积分值，当然这种估计(或近似)是非常粗略的。
     按照图中的提示，求出上述面积的数学期望，就完成了蒙特卡洛积分。
4.蒙特卡洛方法python实例
首先看一个经典的用蒙特卡洛方法求π" role="presentation" style=" box-sizing: border-box; outline: 0px; font-family: "Microsoft YaHei", "SF Pro Display", Roboto, Noto, Arial, "PingFang SC", sans-serif; display: inline; line-height: normal; font-size: 16px; text-align: left; overflow-wrap: break-word; float: none; direction: ltr; max-width: none; max-height: none; min-width: 0px; min-height: 0px; border-width: 0px; border-style: initial; border-color: initial; color: rgb(77, 77, 77); font-variant-ligatures: common-ligatures; background-color: rgb(255, 255, 255);  ">π值。
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
# 投点次数
n = 10000
# 圆的信息
r = 1.0         # 半径
a, b = (0., 0.) # 圆心
# 正方形区域边界
x_min, x_max = a-r, a+r
y_min, y_max = b-r, b+r
# 在正方形区域内随机投点
x = np.random.uniform(x_min, x_max, n) # 均匀分布
y = np.random.uniform(y_min, y_max, n)
# 计算 点到圆心的距离
d = np.sqrt((x-a)**2 + (y-b)**2)
# 统计 落在圆内的点的数目
res = sum(np.where(d # 计算 pi 的近似值(Monte Carlo方法的精髓：用统计值去近似真实值)
pi = 4 * res / n
print('pi: ', pi)
# 画个图看看
fig = plt.figure() 
axes = fig.add_subplot(111) 
axes.plot(x, y,'ro',markersize = 1)
plt.axis('equal') # 防止图像变形
circle = Circle(xy=(a,b), radius=r, alpha=0.5)
axes.add_patch(circle)
plt.show()
简单介绍下思路：   正方形内部有一个相切的圆，它们的面积之比是π/4。现在，在这个正方形内部，随机产生n个点，计算它们与中心点的距离，并且判断是否落在圆的内部。若这些点均匀分布，则圆周率 pi=4 * res/n, 其中res表示落到圆内投点数 n:表示总的投点数。
然后看一个求定积分的例子。
import numpy as np
import matplotlib.pyplot as plt

'''蒙特卡罗方法求函数 y=x^2 在[0,1]内的定积分(值)'''
def f(x):
    return x**2

# 投点次数
n = 10000

# 矩形区域边界
x_min, x_max = 0.0, 1.0
y_min, y_max = 0.0, 1.0     

# 在矩形区域内随机投点
x = np.random.uniform(x_min, x_max, n) # 均匀分布
y = np.random.uniform(y_min, y_max, n)

# 统计 落在函数 y=x^2图像下方的点的数目
res = sum(np.where(y 1, 0))

# 计算 定积分的近似值(Monte Carlo方法的精髓：用统计值去近似真实值)
integral = res / n

print('integral: ', integral)

# 画个图看看
fig = plt.figure() 
axes = fig.add_subplot(111) 
axes.plot(x, y,'ro',markersize = 1)
plt.axis('equal') # 防止图像变形

axes.plot(np.linspace(x_min, x_max, 10), f(np.linspace(x_min, x_max, 10)), 'b-') # 函数图像
#plt.xlim(x_min, x_max)

plt.show()