UA MATH575B 数值分析下计算统计物理例题2

理论解法
- C-K方程法
- 特征值法（近似解）
模拟解法
- Rejection Sampling
- Importance Sampling

一个位于原点x0=0x_0=0x0=0的粒子源，它释放出的例子下一个时刻有2/3的概率的停在原地，有1/3的概率往右移动一个单位；当粒子位于x=1,2,⋯,19x=1,2,\cdots,19x=1,2,⋯,19的位置时，有2/3的概率往左移动一个单位，有1/3的概率往右移动一个单位；x=20x=20x=20是一个吸收状态。求这个吸收状态的吸收率γ\gammaγ。

理论解法

这个粒子的运动可以用一个Markov链表示，记概率转移矩阵为TTT：

%% Construct transition matrix
T = zeros(21,21);
T(1,1) = 2/3; T(21,21) = 1;
for i = 1:21for j = 1:21if j == (i-1) && i<21T(i,j) = 2/3;elseif j == (i+1)T(i,j) = 1/3;endend
end

C-K方程法

用Chapman-Kolmogorov方程的思想，记E(x0→A)E(x_0 \to A)E(x0→A)表示粒子从状态x0x_0x0出发转移到状态子集AAA中所需的时间的期望，则
E(x→A)=∑y∈AT(x,y)+∑y∉AT(x,y)(1+E(y→A))E(x \to A) = \sum_{y \in A} T(x,y) + \sum_{y \notin A} T(x,y) ( 1 + E(y \to A))E(x→A)=y∈A∑T(x,y)+y∈/A∑T(x,y)(1+E(y→A))
第一项表示经过一步转移就到了状态子集AAA中的概率，第二项表示如果第一步没能直接转移到位，记转移到了状态yyy，对应的概率为T(x,y),y∉AT(x,y),y\notin AT(x,y),y∈/A，然后E(y→A)E(y \to A)E(y→A)表示粒子从状态yyy出发转移到状态子集AAA中所需的时间的期望，但此时已经走了一步了所以加1，
∑y∈AT(x,y)+∑y∉AT(x,y)(1+E(y→A))=1+∑y∉AT(x,y)E(y→A)\sum_{y \in A} T(x,y) + \sum_{y \notin A} T(x,y) ( 1 + E(y \to A)) = 1+\sum_{y \notin A} T(x,y) E(y \to A)y∈A∑T(x,y)+y∈/A∑T(x,y)(1+E(y→A))=1+y∈/A∑T(x,y)E(y→A)
记Ex=E(x→{20})E_x = E(x \to \{20\})Ex=E(x→{20})，则
E0=1+23E0+13E1E_0 = 1 + \frac{2}{3} E_0 + \frac{1}{3} E_1E0=1+32E0+31E1
当x=1,⋯,18x=1,\cdots,18x=1,⋯,18时，
Ex=1+23Ex−1+13Ex+1E_x = 1 + \frac{2}{3} E_{x-1} + \frac{1}{3} E_{x+1}Ex=1+32Ex−1+31Ex+1
当x=19x = 19x=19时，
E19=1+23E18E_{19} = 1 + \frac{2}{3} E_{18}E19=1+32E18
由此我们得到了一个有二十个未知量、二十个方程的线性方程组，求解这个线性方程组我们可以得到E0=6291390E_0=6291390E0=6291390。也就是说在原点放出一个粒子，它平均需要移动6291390步才能到达吸收状态，因此吸收状态的吸收率是γ=1/6291390\gamma=1/6291390γ=1/6291390。

特征值法（近似解）

假设初始状态的概率分布是π0\pi_0π0，则nnn步后的状态分布是πn=π0Tn\pi_n = \pi_0 T^nπn=π0Tn。做转移概率矩阵的对角化：T^=VΛV−1\hat{T} = V\Lambda V^{-1}T^=VΛV−1，则πn=π0VΛnV−1\pi_n=\pi_0V \Lambda^n V^{-1}πn=π0VΛnV−1。注意到这个Markov链只有一个吸收状态，所以不论π0\pi_0π0取什么值，当nnn足够大时，πn\pi_nπn都会退化为δ(xn−20)\delta(x_n - 20)δ(xn−20)，并且状态202020对应的特征值为0。因此我们可以用第二大的特征值来近似πn\pi_nπn，
πn≈π0λ2n\pi_n \approx \pi_0 \lambda_2^nπn≈π0λ2n
吸收率就近似为1−λ21-\lambda_21−λ2

模拟解法

Rejection Sampling

这是老师给的答案，大概是用python来写，先import一下random的包然后设置初始值。当粒子没有到状态20时做rejection sampling，输出所有例子最快到达吸收状态的时间。

根据这个code可以画出停时的分布

估计停时的均值为6.2546e+6，吸收率是这个均值的倒数。这个采样最大的问题是慢，因为往右走的概率更小。

Importance Sampling

在rejection sampling中，我们用的是原始的Markov链的转移概率，现在考虑做重要性采样加快模拟的速度。这个Markov链比较有意思，相当于是一个不对称的一维的随机游走，所以要构造一个新的Markov链来采样只需要一个参数ppp（code第二行，手动输入）。接下来用这个新的Markov链做rejection sampling，xxx记录的是新的Markov链的状态转移，compensation记录的是importance sampling的权重，因此每次计数需要加compensation乘以1。
从这个结果可以看出，ppp越大吸收率估计量扭曲得越厉害（Markov链被扭曲了），但速度会加快。要在扭曲和速度之间做tradeoff，可以把这两个Markov链合起来，当状态小于5时就用原来的Markov链，大于5时就用参数为ppp的Markov链。

UA MATH575B 数值分析下计算统计物理例题2相关推荐

UA MATH575B 数值分析下计算统计物理例题1
UA MATH575B 数值分析下计算统计物理例题1 统计物理方法的解析解 Markov链理论解数值解 Monte Carlo模拟. 一道有趣的统计物理的题目.下面这个简单的迷宫中,一只老鼠一开 ...
UA MATH575B 数值分析下VI 统计物理的随机模拟方法2
UA MATH575B 数值分析下VI 统计物理的随机模拟方法2 扩散过程的停时问题理论解模拟解扩散过程的停时问题简单一点,考虑一维的随机扩散问题,假设扩散系数D=12D=\frac{1}{2 ...
UA MATH575B 数值分析下VI 统计物理的随机模拟方法1
UA MATH575B 数值分析下VI 统计物理的随机模拟方法1 模拟SDE的解用蒙特卡罗方法估计SDE解的自相关函数模拟SDE的解以最简单的随机微分方程为例,考虑随机游走 dxt=ξtdtdx ...
UA MATH575B 数值分析下统计物理的随机模拟方法5
UA MATH575B 数值分析下统计物理的随机模拟方法5 Ising Model Gibbs Sampling Glauber Dynamics 这一讲介绍Ising Model,它是MCMC与G ...
UA MATH575B 数值分析下统计物理的随机模拟方法4
UA MATH575B 数值分析下统计物理的随机模拟方法4 这一讲介绍MCMC方法,这个方法最早出现在Metropolis在1953年发在J Chem Phys上的Equation of state ...
UA MATH575B 数值分析下IV 带约束的优化
UA MATH575B 数值分析下IV 带约束的优化问题带等式约束的优化问题带不等式约束的优化问题同时带等式约束与不等式约束的优化问题今天不想敲公式,就不写理论了,反正方法也就是前面的Newt ...
UA MATH575B 数值分析下III 图像恢复
UA MATH575B 数值分析下III 图像恢复图像去噪的优化模型算法实现图像去噪的优化模型假设VVV代表原图,观测到的图像是WWW,WWW是原图与噪声的叠加W=V+ξW=V+\xiW=V+ ...
UA MATH575B 数值分析下I 梯度下降
UA MATH575B 数值分析下I 梯度下降梯度下降最速下降下降算法实现举例梯度下降的应用对于凸优化问题min⁡xf(x)\min_x f(x)minxf(x),最主流的数值计算方法是下 ...
UA MATH575B 数值分析下II 牛顿算法
UA MATH575B 数值分析下II 牛顿算法 Pure Newton算法 Damped Newton算法 Levenberg-Marquardt算法 Quasi-Newton算法割线法 Broy ...

UA MATH575B 数值分析下计算统计物理例题2

UA MATH575B 数值分析下计算统计物理例题2

理论解法

C-K方程法

特征值法（近似解）

模拟解法

Rejection Sampling

Importance Sampling

UA MATH575B 数值分析下计算统计物理例题2相关推荐

最新文章

热门文章

UA MATH575B 数值分析下 计算统计物理例题2

UA MATH575B 数值分析下 计算统计物理例题2

理论解法

C-K方程法

特征值法（近似解）

模拟解法

Rejection Sampling

Importance Sampling

UA MATH575B 数值分析下 计算统计物理例题2相关推荐

最新文章

热门文章

UA MATH575B 数值分析下计算统计物理例题2

UA MATH575B 数值分析下计算统计物理例题2

UA MATH575B 数值分析下计算统计物理例题2相关推荐