0

我正在使用一个简单的二元正态模型，该模型具有一些非常规的先验。我遇到的主要问题是我的后验结果从一次运行到下一次不一致，我猜这与连续样本之间的高度相关性问题有关。这是我的具体问题。

  1. 获得 N 个独立样本的最佳方法是什么？目前，我一直在调用 sample() 来获得一个大链（例如长度 10,000），然后从 1,000 开始每 100 个样本。但是现在看看其中一个参数的自相关曲线，看起来我需要至少每 500 个样本进行一次！（我还可以使用互信息来更好地了解滞后之间的依赖关系。）

  2. 我一直在遵循pymc3 教程中随机波动率示例中描述的拟合过程。特别是我首先找到 MAP，然后用它生成一个 NUTS() 对象，然后取一个简短的样本，然后使用它生成另一个 NUTS() 对象，使用 gamma=0.25 (???)，最后得到我的大样本。我不知道这是否合适或者我是否需要 gamma=0.25。

  3. 此外，在同一个示例中，还有指数分布的测试值。我不知道我是否需要这些。（默认使用均值有什么问题？）

这是我正在使用的实际模型。

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])
4

1 回答 1

0

我不确定你在用你建立的复杂的先验结构做什么，但我认为那里有问题。

我将模型简化为：

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)  # I have actual data!

with pymc.Model():
    corr = Uniform('corr', lower=0, upper=1)
    corr_inv = th.stack([th.stack([1, -corr]),
                         th.stack([-corr, 1])]) / (1 - th.sqr(corr))
    mu = Normal('mu', mu=0, sd=1, shape=2)

    MvNormal('data', 
             mu=mu,
             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])
    trace2 = pymc.sample(10000, step=step2, start=trace1[-1])

具有很大的收敛性。我认为您也可以删除 gamma 参数。

于 2015-12-11T12:29:24.600 回答