使用PyMC的最简单线性模型

Question

假设我想用一些数据来估算一个简单的 y= m * x 问题的斜率：

x_data = np.array([0,1,2,3])
y_data = np.array([0,1,2,3])

很明显，斜率是 1。但是，当我在 PyMC 中运行这个时，我得到了 10

slope  = pm.Uniform('slope',  lower=0, upper=20)

@pm.deterministic
def y_gen(value=y_data, x=x_data, slope=slope, observed=True):
  return slope * x

model = pm.Model([slope])
mcmc  = pm.MCMC(model)
mcmc.sample(100000, 5000)

# This returns 10
final_guess = mcmc.trace('slope')[:].mean()

但它应该是 1！

注意：以上是使用 PyMC2 的情况。

Answer 1

你需要为可能性设置 value=y_data, observed=True。另外，有一点小建议，你不需要创建一个模型对象。只需把你的节点（或者调用 locals()）传给 MCMC 就可以了。

Answer 2

你需要定义一个可能性，试试这个：

import pymc as pm
import numpy as np

x_data = np.linspace(0,1,100)
y_data = np.linspace(0,1,100)

slope  = pm.Normal('slope',  mu=0, tau=10**-2)
tau    = pm.Uniform('tau', lower=0, upper=20)

@pm.deterministic
def y_gen(x=x_data, slope=slope):
  return slope * x

like = pm.Normal('likelihood', mu=y_gen, tau=tau, observed=True, value=y_data)

model = pm.Model([slope, y_gen, like, tau])
mcmc  = pm.MCMC(model)
mcmc.sample(100000, 5000)

# This returns 10
final_guess = mcmc.trace('slope')[:].mean()

它返回10，因为你只是从一个均匀的先验分布中抽样，而10就是这个分布的期望值。