问题描述
完整代码可在此 PyMC2笔记本中获得,并在此 PyMC3笔记本中.
The full code is available in this notebook for PyMC2 and in this notebook for PyMC3.
PyMC2的相关代码是
The relevant code for PyMC2 is
def analyze(data):
# priors might be adapted here to be less flat
mu = pymc.Normal('mu', 0, 0.000001, size=2)
sigma = pymc.Uniform('sigma', 0, 1000, size=2)
rho = pymc.Uniform('r', -1, 1)
@pymc.deterministic
def precision(sigma=sigma,rho=rho):
ss1 = float(sigma[0] * sigma[0])
ss2 = float(sigma[1] * sigma[1])
rss = float(rho * sigma[0] * sigma[1])
return np.linalg.inv(np.mat([[ss1, rss], [rss, ss2]]))
mult_n = pymc.MvNormal('mult_n', mu=mu, tau=precision, value=data.T, observed=True)
model = pymc.MCMC(locals())
model.sample(50000,25000)
我上面的代码到PyMC3的端口如下:
My port of the above code to PyMC3 is as follows:
def precision(sigma, rho):
C = T.alloc(rho, 2, 2)
C = T.fill_diagonal(C, 1.)
S = T.diag(sigma)
return T.nlinalg.matrix_inverse(T.nlinalg.matrix_dot(S, C, S))
def analyze(data):
with pm.Model() as model:
# priors might be adapted here to be less flat
mu = pm.Normal('mu', mu=0., sd=0.000001, shape=2, testval=np.mean(data, axis=1))
sigma = pm.Uniform('sigma', lower=1e-6, upper=1000., shape=2, testval=np.std(data, axis=1))
rho = pm.Uniform('r', lower=-1., upper=1., testval=0)
prec = pm.Deterministic('prec', precision(sigma, rho))
mult_n = pm.MvNormal('mult_n', mu=mu, tau=prec, observed=data.T)
return model
model = analyze(data)
with model:
trace = pm.sample(50000, tune=25000, step=pm.Metropolis())
PyMC3版本可以运行,但是显然不会返回预期的结果.任何帮助将不胜感激.
The PyMC3 version runs, but clearly does not return the expected result. Any help would be highly appreciated.
推荐答案
pymc.Normal的呼叫签名为
The call signature of pymc.Normal is
In [125]: pymc.Normal?
Init signature: pymc.Normal(self, *args, **kwds)
Docstring:
N = Normal(name, mu, tau, value=None, observed=False, size=1, trace=True, rseed=True, doc=None, verbose=-1, debug=False)
请注意,pymc.Normal
的第三个位置参数是tau
,而不是标准偏差sd
.
Notice that the third positional argument of pymc.Normal
is tau
, not the standard deviation, sd
.
因此,由于pymc
代码使用
mu = Normal('mu', 0, 0.000001, size=2)
应使用相应的pymc3
代码
mu = pm.Normal('mu', mu=0., tau=0.000001, shape=2, ...)
或
mu = pm.Normal('mu', mu=0., sd=math.sqrt(1/0.000001), shape=2, ...)
自tau = 1/sigma**2
起.
通过这一更改,您的pymc3代码会生成(类似)
With this one change, your pymc3 code produces (something like)
这篇关于与PyMC3的贝叶斯相关的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持!