Monte Carlo

Orbifold

Last updated on Jan 24, 2024 4 min read

随机数生成

线性同余法

若要生成 $[0, 1]$ 上的随机数，可通过如下递推式生成伪随机数：

$r_{i + 1} = (a r_{i} + c) \mod M$

所得 $r_{i} / M$ 即想要的随机数。

a=5;c=1;r1=3;m=2^8
N=1000;%随机数个数
r=zeros(1,N);
for i=1:N
    r(i)=mod(a*r1+c,m);
    r1=r(i);
end
figure; plot(1:N,r/m,'-*')

MC 积分

直接抽样：考虑 $x_{i} \sim rand (0, 1)$ 服从随机分布，积分等价于求 $f (x_{i})$ 的期望 $I = \int_{0}^{1} f (x) d x = lim_{N \to \infty} \frac{1}{N} \sum_{i = 1}^{N} f (x_{i})$

缺点：往往需要较大的样本量才能达到较小的误差。

重要抽样

考虑变量 $x_{i} \sim w (x)$ ,其中： $\int_{0}^{1} w (x) d x = 1$

利用：

$y (x) = \int_{0}^{x} w (x^{'}) d x^{'}$

则原积分可以表示为：

$I = \int_{0}^{1} \frac{f (x (y))}{w (x (y))} d y = lim_{N \to \infty} \frac{1}{N} \sum_{i = 1}^{N} \frac{f (x (y_{i}))}{w (x (y_{i}))}$

其中 $y_{i} \sim rand (0, 1)$ ，如果 $\frac{f}{w} \sim const$ 则抽样良好。

难点：构造合适的分布 $w (x)$ 与求反函数 $x (y)$ 。

特定随机分布数构造

Gauss 分布

$w (x) = \frac{1}{\sqrt{2 π}} e^{\frac{1}{2} x^{2}}$

由概率论：

$\frac{1}{12} \overset{―}{rand (0, 1)} - 6 \sim N (0, 1)$

%Gauss分布的随机数；

N=100000;
Gauss=zeros(1,N);
for i=1:12
    Gauss=Gauss+rand(1,N);
end

Gauss=Gauss-6;
figure;plot(Gauss,1:N,'*')

%------判断分布是否为Gauss分布-----------------
Num=501;                                                %区间大小
Max=max(abs(Gauss))+10^-9;                       %边界
y=zeros(1,Num-1);
x=linspace(-Max,Max,Num);   %-Max,Max间分为501个点，500段
h=x(2)-x(1);
for i=1:length(Gauss)
    index=fix((Gauss(i)+Max)/h)+1;
    y(index)=y(index)+1;
end
figure;bar(x(1:end-1),y);                       %生成数绘图

Von Neumanm舍选法

考虑 $[a, b]$ 区间上的有界分布 $w (x)$ , $w (x) \leq c \equiv w^{'} (x)$ ,则可如下生成按 $w (x)$ 分布的随机数：

生成 $[a, b]$ 内的均匀分布随机数 $x$ $x = (b - a) rand (0, 1) + a$
生成 $[0, c]$ 内的均匀分布随机数 $y$ $y = c rand (0, 1)$
若 $y \leq w (x)$ ，接受随机数 $x$ w为所需 $w (x)$ 分布的随机数，否则重新抽取 $(x, y)$

%利用Von Neumann舍选法产生随机数
a=0;b=1.0;
h=0.001;
w=@(x) 6/5*(1-0.5*x.^2);%随意选取的分布
N=5*10^6;
xx_s=Von_Neumann(w,a,b,h,N);
plot(xx_s(1:N),1:N,'*')

%---------统计随机数的分布------------
M=1000;
yy1=zeros(1,M);
h1=(b-a)/M;
xx1=a:h1:b;
for i=1:length(xx_s)
    index=fix(xx_s(i)/h1)+1;
    yy1(1,index)=yy1(1,index)+1;
end
figure;bar(xx1(1:end-1),yy1);                       %生成数绘图
S=interg(a,b,w);
yy1=yy1/N/h1*S; %因为yy1/N/h1在[a,b]上的积分为1，为了与w比较，乘以其在[a,b]上的积分S
figure;plot(xx1(1:M),yy1,'r*-');title('Von Neumann舍选法产生的随机数分布');
hold on
plot(xx1(1:M),w(xx1(1:M)),'b')
legend('Von Neumann舍选法产生的随机数分布','精确分布')

%生成步长为$h$的[a,b]区间N个按w(x)分布的随机数函数
function xx_s=Von_Neumann(w,a,b,h,N)
if nargin < 5        %nargin是用来判断输入变量个数的函数 
    h = 10^-4; %步长标准
end
x=a:h:b;
y_m=max(w(x))*1.001;%选择的w'
i=0;     
n=0;
xx_s=zeros(N,1);
while i<N
    xx=(b-a)*rand(1,1)+a;%产生[a,b]区间上的均匀分布的随机数
    yy=y_m*rand(1,1);
    n=n+1;   %统计产生了多少次随机数
    if yy<=w(xx) %舍选法的判断标准
        i=i+1;
        xx_s(i)=xx; 
    end
end
end

function S=interg(a,b,w) %w函数在[a,b]上的积分
    n=10000;
    h=(b-a)/n;
    xx_h=a:h:b;
    ww=w(xx_h);
    S=(4*sum(ww(2:2:n))+2*sum(ww(3:2:n-1))+ww(1)+ww(n+1))*h/3; %simpson积分公式
end

Metropolis算法

Metropolis算法与辐射原理一致，设变量序列 $X \sim w (X)$ ，应决定 $X_{n}$ 的下一点 $X_{n + 1}$ 点生成。

将 $w (X)$ 类比为能量，若 $X_{t}$ 相对于 $X_{n}$ 处于高能级, 即 $w (X_{t}) > w (X_{n})$ ，则 $X_{t} \to X_{n}$ 的自发辐射概率为1, $X_{t}$ 总是可接受的；
若 $X_{t}$ 相对于 $X_{n}$ 处于低能级, 即 $w (X_{t}) < w (X_{n})$ ，则 $X_{t} \to X_{n}$ 的受激辐射概率为 $r = \frac{w (X_{t})}{w (X_{n})}$ ,仅有概率 $r$ 使得 $X_{t}$ 是可接受的；

一般严格算法如下：

取 t=0对应无迭代初始状态;
利用先验分布 $π^{(0)}$ 生成初始元素 $x_{0}$ ;
进行如下循环直到 $t = M$ ， $M$ 为采样数；
- $t = t + 1$ ;
- 利用建议分布 $q (x ∣ x_{t - 1})$ 生成建议点 $x *$ ；
- 计算接受概率 $α = \min (1, \frac{p (x^{*})}{p (x^{(t - 1)})})$
- 生成均匀随机数 $u \sim rand (0, 1)$
- 若 $u \leq α$ ，则接受建议点 $x_{t} = x *$ ;
- 否则 $x_{t} = x_{t - 1}$

通常取 $π^{(0)} \sim N (0, 1)$ ，而 $q (x_{t} ∣ x_{t - 1}) \sim N (x_{t - 1}, 1)$ .


clear; clc;
close all;

N=50000;
a=0;
b=10;
w = @(x) (1 + x.^2).^-1;
xx=mymetro(w,N);
figure;
histogram(xx,200,'Normalization','pdf');%绘制 x 的概率密度函数估算值。


Num=501;                                                %区间大小
Max=max(abs(xx))+10^-9;                       %边界
y=zeros(1,Num-1);
x=linspace(-Max,Max,Num);
h=x(2)-x(1);
for i=1:length(xx)
    index=fix((xx(i)+Max)/h)+1; %判定xx(i)处于哪个子区间，fix为舍去小数取整运算
    y(index)=y(index)+1;          %该子区间的随机数数目+1
end
figure;bar(x(1:end-1),y); 


%精确值绘图
figure,fplot(w,[-Max,Max],'r');
hold on;
S=interg(-Max,Max,w);
y=y/length(xx)/h*S;  %生成的随机数的数目为n(length(xx))，走的总步数为N
plot(x(1:Num-1),y,'*');title('验证随机数分布');
legend('精确分布','随机数分布')
hold off


function S=interg(a,b,w) %w函数在[a,b]上的积分 
    n=10000;
    h=(b-a)/n;
    xx_h=a:h:b;
    ww=w(xx_h);
    %S=sum(w(xx_h).*(b-a)/n)-w(a)*(b-a)/2/n-w(b)*(b-a)/2/n; %梯形算法
    S=(4*sum(ww(2:2:n))+2*sum(ww(3:2:n-1))+ww(1)+ww(n+1))*h/3; %simpson积分公式
end

function x=mymetro(w,N)
    x=zeros(1,N);
    x(1)=randn;%标准正态分布生成初始随机数
    t = 1; % 抽样编号
    while t<N
        t=t+1;
        %建议采样
        xs= normrnd(x(t-1) ,1); % 按照建议抽样函数进行采样
        r=min([1, w(xs)/w(x(t-1))]); % 接受概率
        %判断是否接受
        u = rand; % Metropolis 采样的核心部分
        if u < r
            x(t) = xs;
        else
            x(t) = x(t-1);
        end  
    end
end

计算物理

Monte Carlo

随机数生成

MC 积分

特定随机分布数构造

Orbifold

Undergraduated Student