在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$上时，包括雅可比。

大都会黑斯廷斯接受-拒绝步骤不正确！

if 0 < probab:

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