序列蒙特卡洛建模 python

主要解决仿真模拟的问题。

什么是蒙特卡罗方法。

• 蒙特卡罗方法又称为统计模拟法、随机抽样技术，是一种随机模拟方法，以概率和统计理论方法为基础的一种计算方法，是使用随机数或者伪随机数来解决很多计算问题的方法。将所求解的问题同一定的概率模型相联系，用电子计算机实现统计模拟或抽样，以获得问题的近似解。

基本思想

• 当所求的问题的解是某个事件的概率，或者是某个随机变量的数学期望，或者是与概率，数学期望有关的量时，通过某种试验的方法，得到该事件发生的概率,或者该随机变量若干个具体观察值的算术平方根，通过它得到问题的解。
• 当随机变量的取值仅为1或者0时，他的数学期望就是某个事件的概率。或者说某种事件的概率也是随机变量（仅取值1或者0）的数学期望。

例如上图所示，可以通过大量的随机模拟投掷来得到每个区域的概率。这样就可以转化到随机边界或者不定边界的问题。例如下图，可以求解具体A、B两个图形的面积。将求解面积的问题转化成求解概率的问题。

蒙特卡洛方法的特点

优点：

• 能够比较逼真地描述具有随机性质的事物的特点及物理实验过程。
• 受几何条件限制小
• 收敛速度与问题的维数无关
• 具有同时计算多个方案与多个未知量的能力
• 误差容易确定
• 程序结构简单，易于实现。
  缺点：
• 收敛速度慢，耗时比较多
• 误差具有概率性
• 在粒子输运问题中，计算结果与系统大小有关。

所以在使用蒙特卡洛方法的时候，要扬长避短，只对问题中难以使用解析（或者数值）方法处理的部分使用蒙塔卡罗方法，对那些能使用解析解处理的方法要尽量使用解析方法

主要应用范围：

• 粒子输运问题（实验物理，反应堆物理）
• 统计物理
• 典型数学问题
• 真空技术
• 激光技术
• 医学
• 生物
• 探矿等

什么是随机数

• 在连续型随机变量的分布中，最简单而且最基本的分布是单位均匀分布。由该分布抽取的简单子样称为随机数序列，其中每一个体称为随机数
• 符号：
• 两个特点：独立性和均匀性
  产生随机数：
  随机数表：
• 随机数表是由0123456789十个数字组成，每个数字以0.1的等概率出现，数字之间相互独立，这些数字序列叫做随机数序列。
• 如果要得到n位有效数字的随机数。只需要将表中每n个相邻的随机数字合并在一起，且在最高位的前边加上小数点即可。例如，某随机数表的第一行数字为7 6 3 4 2 5 8 9 1 ，要想得到三位有效数字的随机数一次为0.763 0.425 0.891 .。。。如果是四位：则0.7634 0.2589.。。。

物理方法：

• 利用某些物理现象，在计算机上增加些特殊设备，可以在计算机上直接产生随机数。
• 作为随机数发生器的物理源主要有两种：一种是根据放射性物质的放射性，另一种是利用计算机的固有噪声。
• 一般情况下，任意一个随机数在计算机内总是用二进制的数表示的：
  
  其中 或者为0，或者为1。因此，利用物理方法在计算机产生随机数，就是要产生只取0或1的随机数字序列，数字之间相互独立，每个数字取0或1的概率均为0.5
  缺点：
• 随机数表需在计算机中占有很大内存，而且也难以满足蒙特卡罗方法对随机数需要量非常大的要求，因此，该方法不适于在计算机上使用。
• 物理方法产生的随机数序列无法重复实现，不能进行程序复算。给验证结果带来很大困难。而且增加随机数发生器和电路联接等附加设备，费用昂贵。因此该方法也不适合在计算机上使用。

伪随机数

用递推公式

产生随机数列：

伪随机数存在的两个问题：

• 递推公式和初始值确定后，整个随机数序列便被唯一确定。不满足随机数相互独立的要求。
• 由于随机数序列是由递推公式确定的，而在计算机上所能表示的【0，1】上的数又是有限的，因此，这种方法产生的随机数序列就不可能不出现重复。随机数序列出现周期性的循环现象。
  解决方法：
• 不能从本质上加以改变，但只要递推公式选的比较好，随机数间的相互独立性是可以近似满足的。
• 因为使用蒙塔卡罗方法解任何具体问题时，所使用的随机数的个数总是有限的，只要所使用的随机数的个数不超过伪随机数序列出现循环现象时的长度就可以了。

应用：蒙塔卡罗方法计算积分

可以通俗地说，蒙特卡罗方法是用随机试验的方法计算积分，即将所要计算的积分看作服从某种分布密度函数f®的随机变量ｇ®的数学期望 。

通过某种试验，得到Ｎ观察值r1，r2，…，rN（用概率语言来说，从分布密度函数f®中抽取Ｎ个子样r1，r2，…，rN，），将相应的Ｎ个随机变量的值g(r1)，g(r2)，…，g(rN)的算术平均值

例：求0-b的积分（求面积）：

那么这个积分是否可以用一个近似的长方形来表示呢？是可以的

根据中值定理一定有一个点b*f(x)=那个积分。

积分问题可以抽象出一个高，这个高就是那个平均值。

蒙特卡洛方法步骤如下：

• 1在区间[a,b]上利用计算机均匀产生n个随机数x1,x2…xn，使用matlab软件的unifrnd命令实现。
• 2计算每一个随机数想对应的被积函数值f(x1),f(x2),f(xn)
• 计算被积函数值的平均值。



• 4原式约等于

蒙特卡罗方法：随机投点试验求近似解

例子：

给定曲线y =2 – x2 和曲线y3 = x2，曲线的交点为：P1( – 1，1 )、P2( 1，1 )。曲线围成平面有限区域，用蒙特卡罗方法计算区域面积。

P=rand(10000,2);
x=2*P(:,1)-1;
y=2*P(:,2);
II=find(y<=2-x.^2&y.^3>=x.^2);
M=length(II);
S=4*M/10000
plot(x(II),y(II),'g.')

例5.14 计算

其中D为y= x – 2与y2 = x 所围D的边界曲线交点为：(–1，1)，(4，2)，被积函数在求积区域内的最大值为16。积分值是三维体积，该三维图形位于立方体区域

0≤ x ≤4，–1≤ y ≤2，0 ≤ z ≤16

内，立方体区域的体积为192。

data=rand(10000,3);
 x=4*data(:,1);
 y=-1+3*data(:,2);
 z=16*data(:,3);
 II=find(x>=y.^2&x<=y+2&z<=x.*(y.^2));
 M=length(II);
 V=192*M/10000

P=rand(10000,4);
 x=-1+2*P(:,1); y=-1+2*P(:,2);
 z=P(:,3);u=2*P(:,4);
 II=find(z>sqrt(x.^2+y.^2)&z<=1&u<=x.^2+y.^2+z.^2);
 M=length(II);
 V=8*M/10000

function icecream(m,n) 
if nargin==0,m=20;n=100;end
t=linspace(0,2*pi,n);
r=linspace(0,1,m);
x=r'*cos(t);y=r'*sin(t);
z1=sqrt(x.^2+y.^2);
z2=1+sqrt(1+eps-x.^2-y.^2);
X=[x;x];Y=[y;y];
Z=[z1;z2];
mesh(X,Y,Z)
view(0,-18)
colormap([0 0 1]),axis off

x1=-1:0.1:1; 
y1=x1.^2.^(1/3); 
x2=1:-0.1:-1;
 y2=2-x2.^2; 
fill([x1,x2],[y1,y2],'g')
axis off