1

我不确定这算作问题还是错误报告。我在这里发布了一个 GitHub 要点：https ://gist.github.com/jbwhit/a9012e04b0f48e582c22

我发现这个问题（pymc3：具有多个观察变量的层次模型）是我自己的层次模型的一个很好的起点，但是当我试图以任何实质性的方式修改它时遇到了困难。

首先，有效的模型和设置：

import numpy as np
import pymc3 as pm

n_individuals = 200
points_per_individual = 10
means = np.random.normal(30, 12, n_individuals)
observed = np.random.normal(means, 1, (points_per_individual, n_individuals))

model = pm.Model()
with model:
    hyper_mean = pm.Normal('hyper_mean', mu=0, sd=100)
    hyper_sigma = pm.HalfNormal('hyper_sigma', sd=3)

    means = pm.Normal('means', mu=hyper_mean, sd=hyper_sigma, shape=n_individuals)
    sigmas = pm.HalfNormal('sigmas', sd=100)

    ye = pm.Normal('ye', mu=means, sd=sigmas, observed=observed)

    trace = pm.sample(10000)

以上所有工作都按预期工作（并且痕迹看起来不错）。下一段代码进行了一次更改（将 T 分布替换为正态分布）：

model = pm.Model()
with model:
    hyper_mean = pm.Normal('hyper_mean', mu=0, sd=100)
    hyper_sigma = pm.HalfNormal('hyper_sigma', sd=3)

    ### Changed to a T distribution ###
    means = pm.StudentT('means', nu=hyper_mean, sd=hyper_sigma, shape=n_individuals)

    sigmas = pm.HalfNormal('sigmas', sd=100)

    ye = pm.Normal('ye', mu=means, sd=sigmas, observed=observed)

    trace = pm.sample(10000)

以下是输出：

Assigned NUTS to hyper_mean
Assigned NUTS to hyper_sigma_log
Assigned NUTS to means
Assigned NUTS to sigmas_log
---------------------------------------------------------------------------
PositiveDefiniteError                     Traceback (most recent call last)
<ipython-input-12-69f59e2f3d47> in <module>()
     18     ye = pm.Normal('ye', mu=means, sd=sigmas, observed=observed)
     19 
---> 20     trace = pm.sample(10000)

/Users/jonathan/miniconda2/envs/pymc3/lib/python3.5/site-packages/pymc3/sampling.py in sample(draws, step, start, trace, chain, njobs, tune, progressbar, model, random_seed)
    121     """
    122     model = modelcontext(model)
--> 123 
    124     step = assign_step_methods(model, step)
    125 

/Users/jonathan/miniconda2/envs/pymc3/lib/python3.5/site-packages/pymc3/sampling.py in assign_step_methods(model, step, methods)
     66             selected_steps[selected].append(var)
     67 
---> 68     # Instantiate all selected step methods
     69     steps += [s(vars=selected_steps[s]) for s in selected_steps if selected_steps[s]]
     70 

/Users/jonathan/miniconda2/envs/pymc3/lib/python3.5/site-packages/pymc3/sampling.py in <listcomp>(.0)
     66             selected_steps[selected].append(var)
     67 
---> 68     # Instantiate all selected step methods
     69     steps += [s(vars=selected_steps[s]) for s in selected_steps if selected_steps[s]]
     70 

/Users/jonathan/miniconda2/envs/pymc3/lib/python3.5/site-packages/pymc3/step_methods/nuts.py in __init__(self, vars, scaling, step_scale, is_cov, state, Emax, target_accept, gamma, k, t0, model, profile, **kwargs)
     76 
     77 
---> 78         self.potential = quad_potential(scaling, is_cov, as_cov=False)
     79 
     80         if state is None:

/Users/jonathan/miniconda2/envs/pymc3/lib/python3.5/site-packages/pymc3/step_methods/quadpotential.py in quad_potential(C, is_cov, as_cov)
     33         return QuadPotential_SparseInv(C)
     34 
---> 35     partial_check_positive_definite(C)
     36     if C.ndim == 1:
     37         if is_cov != as_cov:

/Users/jonathan/miniconda2/envs/pymc3/lib/python3.5/site-packages/pymc3/step_methods/quadpotential.py in partial_check_positive_definite(C)
     56     if len(i):
     57         raise PositiveDefiniteError(
---> 58             "Simple check failed. Diagonal contains negatives", i)
     59 
     60 

PositiveDefiniteError: Scaling is not positive definite. Simple check failed. Diagonal contains negatives. Check indexes [202]

关于如何让它工作的任何建议？

4

2 回答 2

4

正如我在评论中提到的，尝试运行：

model = pm.Model()
with model:
    hyper_mean = pm.Normal('hyper_mean', mu = 0, sd = 100)
    hyper_sigma = pm.HalfNormal('hyper_sigma', sd = 3)
    nu = pm.Exponential('nu', 1./10, testval = 5.)

    ### Changed to a T distribution ###
    means = pm.StudentT('means', nu = nu, mu = hyper_mean, sd = hyper_sigma, shape = n_individuals)

    sigmas = pm.HalfNormal('sigmas', sd = 100)

    ye = pm.Normal('ye', mu = means, sd = sigmas, observed = observed)

    trace = pm.sample(10000)

换句话说：使用该方法的mu论点pm.StudentTforhyper_mean和nufor 自由度。

一旦它开始工作，您也可以尝试添加该pm.find_MAP方法（如@Chris Fonnesbeck 建议的那样）。

于 2016-07-21T08:17:09.630 回答
4

尝试找到 MAP 估计值并将其用作 MCMC 运行的起点：

start = pm.find_MAP()
trace = pm.sample(10000, start=start)
于 2016-03-04T15:35:42.747 回答