我正在使用一个简单的二元正态模型,该模型具有一些非常规的先验。我遇到的主要问题是我的后验结果从一次运行到下一次不一致,我猜这与连续样本之间的高度相关性问题有关。这是我的具体问题。
获得 N 个独立样本的最佳方法是什么?目前,我一直在调用 sample() 来获得一个大链(例如长度 10,000),然后从 1,000 开始每 100 个样本。但是现在看看其中一个参数的自相关曲线,看起来我需要至少每 500 个样本进行一次!(我还可以使用互信息来更好地了解滞后之间的依赖关系。)
我一直在遵循pymc3 教程中随机波动率示例中描述的拟合过程。特别是我首先找到 MAP,然后用它生成一个 NUTS() 对象,然后取一个简短的样本,然后使用它生成另一个 NUTS() 对象,使用 gamma=0.25 (???),最后得到我的大样本。我不知道这是否合适或者我是否需要 gamma=0.25。
此外,在同一个示例中,还有指数分布的测试值。我不知道我是否需要这些。(默认使用均值有什么问题?)
这是我正在使用的实际模型。
import pymc3 as pymc
import numpy as np
import theano.tensor as th
from pymc3.distributions.continuous import Gamma, Uniform, Normal, Bounded
from pymc3.distributions.multivariate import MvNormal
from pymc3.model import Deterministic
data = np.random.randn(3000, 2) / 300 # I have actual data!
with pymc.Model():
tau = Gamma('tau', alpha=2, beta=1 / 20000)
sigma = Deterministic('sigma', 1 / th.sqrt(tau))
corr = Uniform('corr', lower=0, upper=1)
alpha_sig = Deterministic('alpha_sig', sigma / 50)
alpha_post = Normal('alpha_post', mu=0, sd=alpha_sig)
alpha_pre = Bounded(
'alpha_pre', Normal, alpha_post, np.Inf, mu=0, sd=alpha_sig)
corr_inv = th.stack([th.stack([1, -corr]),
th.stack([-corr, 1])]) / (1 - th.sqr(corr))
MvNormal(
'data', mu=th.stack([alpha_post, alpha_pre]),
tau=tau * corr_inv, observed=data)
map_ = pymc.find_MAP()
step1 = pymc.NUTS(scaling=map_)
trace1 = pymc.sample(1000, step=step1)
step2 = pymc.NUTS(scaling=trace1[-1], gamma=0.25)
trace2 = pymc.sample(10000, step=step2, start=trace1[-1])