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

我正在尝试使用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使用log：

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][:]

最佳答案

当提案为正态分布且支持$（-\ infty，+ \ infty）$时，得到负的标准偏差$ \ sigma $一点也不奇怪。

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

按照书面规定，Metropolis-Hastings接受-拒绝步骤不正确！

if 0 < probab:

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

关于python - 在Metropolis-Hasting中使标准偏差始终为正，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/42225654/