python - 在 Metropolis-Hastings 中使标准差始终为正

Question

我正在尝试使用 Metropolis-Hastings 进行贝叶斯回归。测试数据生成如下(Python代码，我没有复制整个代码):

trueA = 5 ; trueB = 7 ;trueSD = 10 ; sample_size = 261
x = np.arange(-sample_size/8, sample_size/8, (sample_size*2/8)/sample_size)
y = trueA *x + trueB + npr.normal(loc=0, scale=trueSD, size=sample_size)

我将对数似然定义如下:

def likelihood(param):
    a = param[0][0] ; b = param[0][1] ; sd = param[0][2] ; pred = a*x + b
    sumSqError = np.power((y - pred), 2).sum()
    likelihoodsum = ((sample_size/2)*(np.log(1)-np.log(np.power(sd,2)))) + (- 1/(2*np.power(sd,2)) * sumSqError)
    return likelihoodsum

为了说明下一点，我准备了以下函数:

def next_param(param, param_index):
    a_next = param[0][0] ; b_next = param[0][1] ; sd_next = param[0][2]

    if param_index == 0:
        a_next = param[0][0] + npr.normal(0, 0.1)
    elif param_index == 1:
        b_next = param[0][1] + npr.normal(0, 0.1)
    elif param_index == 2:
        sd_next = param[0][2] + npr.normal(0, 0.1)

    return np.array([[a_next, b_next, sd_next]])

这段代码运行良好(接受率足够高，我可以估计参数)，尽管我知道 sd_next 在上面的代码中可能会变成负值，这很奇怪。

因此，我决定对 sd_next 使用日志:

elif param_index == 2:
  sd_next = np.log(param[0][2]) + npr.normal(0, 0.1)

return np.array([[a_next, b_next, np.exp(sd_next)]])

然而，估计的参数与真实值相差甚远。如何使 Metropolis-Hastings 中的标准差始终为正？

JFI，这是 MCMC 部分:

num_sampling = 1000
chain = np.zeros((num_sampling, 1, 3))
chain[0][0][0] = 20 # starting value for a
chain[0][0][1] = 15 # starting value for b
chain[0][0][2] = 15 # starting value for sd

num_accepted = 0
for i in range(num_sampling-1):
    chain_previous = chain[i][:]
    chain_new = np.zeros((1, 1, 3))

    for p in range(3):
        proposal = next_param(chain_previous, p)

        probab = likelihood(proposal) - likelihood(chain_previous)
        if 0 < probab:
            chain_new[0][0][p] = proposal[0][p]
            num_accepted += 1
        else:
            chain_new[0][0][p] = chain[i][0][p]

    chain[i+1] = chain_new[0][:]

Answer 1

当您的提案是正态分布且支持 $(-\infty,+\infty)$ 时，您得到负标准差 $\sigma$ 一点也不奇怪。

Metropolis-Hastings 接受-拒绝步骤还应包括三个参数的先验分布。当提案位于 $\log\sigma$ 时，包括雅可比行列式。

如所写，Metropolis-Hastings 接受-拒绝步骤不正确!

if 0 < probab:

不是接受向建议值移动的正确条件:应该将(对数)概率与(对数)均匀值进行比较。在当前格式中，您会收敛到最大可能性。