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

一道有趣的统计物理的题目。下面这个简单的迷宫中，一只老鼠一开始在3的位置，它在岔路口走每条路的概率是均等的，比如这个时刻在3的位置，那么下一个时刻有1/2的可能在2的位置，有1/2的时刻在7的位置，目标是计算老鼠成功逃出的概率。

统计物理方法的解析解

首先我们可以试图找一下理论解。不妨假设3的位置是一块榴莲，13的位置是一个效果非常好的活性炭除臭盒，12的位置会让气味集中，我们想知道榴莲的味道均匀扩散后（达到平衡状态）迷宫里面臭味的分布情况。根据平衡状态的关系，假设用$r_i$ri​表示平衡状态每个位置的臭味，以位置0为例，
$r_0 = (r_1 + r_4)/2$r0​=(r1​+r4​)/2
除位置12、13外一共有14个这样的方程（因为12为1，13为0），可以用mathematica来解（默默吐槽一下，这个就是算平稳分布而已）：

这些值，比如$r_3=9/23$r3​=9/23，说明位置3的气味分子会有$9/23$9/23流窜到位置12。也就是说小鼠从位置3出发逃出去的概率是9/23=0.391304348…。

Markov链

写出上面那个迷宫的Markov转移概率矩阵：

理论解

假设初始状态的概率分布是$\pi_0$π0​，则$n$n步后的状态分布是$\pi_n = \pi_0 \hat{T}^n$πn​=π0​T^n。做转移概率矩阵的对角化：$\hat{T} = V\Lambda V^{-1}$T^=VΛV−1，则$\pi_n=\pi_0V \Lambda^n V^{-1}$πn​=π0​VΛnV−1。注意到这个Markov链有两个吸收状态（位置12和13），假设平稳分布是$\pi$π，则$\pi = \pi \hat{T}$π=πT^，也就是说吸收状态对应的转移概率矩阵的特征值为1，也就是说$\lim \Lambda^n(13) = 1,\lim \Lambda^n(14) = 1$limΛn(13)=1,limΛn(14)=1。所以
$\pi_{\infty}(13) = (V*V^{-1})(13),\pi_{\infty}(14) = (V*V^{-1})(14)$π∞​(13)=(V∗V−1)(13),π∞​(14)=(V∗V−1)(14)
用python求解

可以得到概率大约是是0.39130435。

数值解

数值计算的思路是用$\pi_n = \pi_0 \hat{T}^n$πn​=π0​T^n做矩阵乘法，当$n$n足够大的时候$\pi_n$πn​会收敛到平稳分布。


这里估计的概率是0.391304348。基于Markov链的这两种方法都还是很准确的。

Monte Carlo模拟.

一个做Monte Carlo模拟的matlab示例。模拟百万次成功逃出的次数是390878，所以估计概率是0.390878。

P = [0 0.5 0 0 0.5 0 0 0 0 0 0 0 0 0 0 0;1/3 0 1/3 0 0 1/3 0 0 0 0 0 0 0 0 0 0; ...
    0 0.5 0 0.5 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0.5 0 0 0 0 0.5 0 0 0 0 0 0 0 0; ...
    1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0.5 0 0 0 0 0 0 0 0.5 0 0 0 0 0 0; ...
    0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0; 0 0 0 1/3 0 0 1/3 0 0 0 0 1/3 0 0 0 0; ...
    0 0 0 0 0 0 0 0 0 0.5 0 0 0.5 0 0 0; 0 0 0 0 0 1/3 0 0 1/3 0 1/3 0 0 0 0 0; ...
    0 0 0 0 0 0 0 0 0 0.5 0 0 0 0 0.5 0; 0 0 0 0 0 0 0 0.5 0 0 0 0 0 0 0 0.5; ...
    0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0; ...
    0 0 0 0 0 0 0 0 0 0 1/3 0 0 1/3 0 1/3; 0 0 0 0 0 0 0 0 0 0 0 0.5 0 0 0.5 0];

rng(0);
Eaten = 0; Escape = 0; 
for i = 1:1000000 % Set free 1 million mice
    x = 3; % Set free 1 million mice at position 3
    while x ~= 12 && x ~= 13 % when the mouse hasn't escaped or been eaten 
        u = rand(1,1);
        sum = P(x+1,1);
        y = 0;
        while sum < u % check which state will y transit to
            sum = sum + P(x+1,y+2);
            y = y + 1;
        end
        x = y;
    end
    switch x
        case 12
            Eaten = Eaten + 1; continue; % record the number of eaten mice
        case 13
            Escape = Escape + 1; continue; % record the number of escaped mice
    end
end