使用 Metropolis 算法的环岛游

机器算法验证 Python 大都会黑斯廷斯
2022-04-08 17:08:08

来自 John Kruschke 的书，第 7 章，Pg。120（为简洁起见）：

一个政治家不断地在一系列岛屿上从一个岛屿到另一个岛屿旅行......他的目标是按比例访问所有岛屿，以便他在人口较多的岛屿上花费更多的时间......他没有甚至确切地知道有多少个岛屿！.. 他的顾问可以向岛上的市长询问他们在岛上有多少人。当政客提议访问相邻的岛屿（向东或向西，通过掷硬币）时，他们可以询问该相邻岛屿的市长该岛上有多少人。

我写了一些 Python 代码，A岛链在哪里run_Metropolis(A, N, start, burn_in)并进行计算。代码被故意简化以使其更易于阅读。

import random
def run_Metropolis(A, N, start, burn_in):

    ''' A: the distribution to be predicted (a list).
        N: the number of trials / moves
        start: starting index in A
        burn_in: the number of  initial moves to be discarded.
    '''

    Posterior = [0] * len(A) # to be returned
    cur_pos = start # starting index

    for i in range(N):

        # proposal distribution is (-1, +1) cordinates, takng into account the list boundaries
        if cur_pos<len(A)-1 and cur_pos>0:
            proposed_positions = [cur_pos-1, cur_pos+1]
        elif cur_pos == len(A)-1:
            #proposed_positions =[cur_pos-1] # if this is enforced, the posterior does not look like A
            proposed_positions = [0,cur_pos-1] 
        elif cur_pos == 0 :
            #proposed_positions = [cur_pos+1] # if this is enforced, the posterior does not look like A
            proposed_positions = [len(A)-1, cur_pos+1]     
        proposed_pos = random.choice(proposed_positions)

        # decide whether you will move
        if A[proposed_pos] > A[cur_pos]: # definitely move
            cur_pos = proposed_pos
        else: # move with a prob proportional to how close A[cur_pos] is to A[proposed_pos] 

            u = random.random() # sample from uniform distribution from 0 to 1

            if float(A[proposed_pos])/A[cur_pos] > u: # move to proposed
                 cur_pos = proposed_pos 
            else: # stay at cur_pos
                pass 

        # if past burn_in, increment the location by 1 in the posterior distribution.
        if i > burn_in: 
            Posterior[cur_pos] +=  1

    return Posterior

 A = [20,30,10,50,80,20,90] # populations of the islands
 posterior_A = run_Metropolis(A, 4000, 3, 1000)
 # Validate the sampling works:
 [n/float(sum(A)) for n in A]
 > [0.0666, 0.1, 0.0333, 0.166, 0.266, 0.066, 0.3] 
 [n/float(sum(posterior_A)) for n in posterior_A]
 > [0.0593, 0.0910, 0.0393, 0.177, 0.263, 0.0666, 0.302]

所以，posterior_A确实反映了人口分布A。Kruschke 没有说清楚的是，为了使这个工作，岛链需要是循环的，即你需要能够从第 7 个岛旅行到第 1 个岛，反之亦然。否则，posterior_A看起来不像A。您可以通过使用注释掉的边界提案分布来测试这一点，这些分布不假设周期性。

问题：

1.这是因为违反了遍历性假设吗？

2.我们能否修改接受标准（移动决策）以解决最东和最西岛屿缺乏对称性的问题？（我猜这是 Metropolis-Hastings）

3.如果可以，能否请您修改代码说明一下？

1个回答

1. 问题不在于遍历性

不，这与遍历性无关。在没有循环的链条中，人们仍然可以从任何一个岛屿移动到任何其他岛屿，并且（假设存在人口差异！）该链条不是周期性的（因为有时会保持不变），因此该链条是遍历的，因此它具有独特的平稳分布。问题在于平稳分布是错误的，即满足详细平衡但不符合预期分布。

2. 大都会-黑斯廷斯

你的想法是正确的。在 Metropolis-Hastings 中，修改接受概率以考虑提案分布中的不对称性：提案被接受的概率为 $\min(r,1)$，其中 \begin{equation} r=\frac{p(x ')}{p(x)}\,\frac{g(x'\rightarrow x)}{g(x \rightarrow x')}。\end{equation} 其中 $p$ 是目标分布，$x$ 是旧状态（岛），$x'$ 是建议状态，$g(y\rightarrow z)$ 是概率（在当链在 $y$ 时提出 $z$ 的连续案例）。在您的循环链中，$g$ 始终是 $1/2$，因此后一个因素消失了，剩下的就是普通的 Metropolis 算法。min(r,1)，其中 其中是目标分布，是旧状态（岛），是建议状态，是当链在。在您的循环链始终为，因此后一个因素消失了，我们剩下的是普通的 Metropolis 算法。

r=p(x′)p(x)g(x′→x)g(x→x′).
pxx′g(y→z)zyg1/2

在新链中，当离开任一边界岛时，$g$ 为 1。在所有其他情况下，$g=0.5$。因此，当旧岛或新岛是边界岛时，必须考虑接受概率中的后一个因素（该因素分别为 0.5 美元或 2 美元）。g当离开任一边界岛时，g为 1。在所有其他情况下g=0.5。因此，当旧岛或新岛是边界岛时，必须考虑接受概率中的后一个因素（该因素分别为0.5或2）。

3. 实施

仅附上相关修改部分（提案和验收）。这会直接计算与 $1$ 不同的特殊情况下 $g$s 的比率。g s 对于特殊情况，它不同于1。

    proposal_ratio = 1 #The ratio of proposal probabilities, change for special cases
        # proposal distribution is (-1, +1) cordinates, takng into account the list boundaries
    if cur_pos<len(A)-1 and cur_pos>0:
        proposed_positions = [cur_pos-1, cur_pos+1]
    elif cur_pos == len(A)-1:
        proposed_positions =[cur_pos-1] 
        proposal_ratio = 0.5
    elif cur_pos == 0 :
        proposed_positions = [cur_pos+1] 
        proposal_ratio = 0.5

    proposed_pos = random.choice(proposed_positions)

    if proposed_pos == 0 or proposed_pos == len(A)-1:
        proposal_ratio = 2

    # decide whether you will move
    r = (float(A[proposed_pos])/A[cur_pos]) * proposal_ratio
    if r>=1: # definitely move
        cur_pos = proposed_pos
    else: # move with prob r 
        u = random.random() # sample from uniform distribution from 0 to 1
        if r > u: # move to proposed
             cur_pos = proposed_pos 
        else: # stay at cur_pos
            pass
其它你可能感兴趣的问题
上一篇预测变量在多元线性回归中的重要性 下一篇MCMC 建模 - 这甚至可以解决吗？