2

我正在将来自 Stan 的 ARM 书中的选定模型“翻译”到 pymc3（我希望尽快将它们发布到 Github 上），并且我有一个关于“不正确的先验”的问题。

我知道 Stan 默认是对参数使用统一的先验。如果无界，这意味着统一（-inf，+inf）

我的问题是：无论如何在pymc3中指定这样的先验？

这是一个示例来说明问题以及我到目前为止所做的尝试。

import numpy as np
import pymc3 as pm

light_speed <-  np.array(28, 26, 33, 24, 34, -44, 27, 16, 40, -2, 29, 22, 24, 21, 25, 
30, 23, 29, 31, 19, 24, 20, 36, 32, 36, 28, 25, 21, 28, 29, 
37, 25, 28, 26, 30, 32, 36, 26, 30, 22, 36, 23, 27, 27, 28, 
27, 31, 27, 26, 33, 26, 32, 32, 24, 39, 28, 24, 25, 32, 25, 
29, 27, 28, 29, 16, 23)

在斯坦（pystan）

 #the model in pystan as specified in stan arm examples

model_string = '''
data {
  int<lower=0> N; 
  vector[N] y;
 }
parameters {
vector[1] beta;
real<lower=0> sigma;
} 
model {
 y ~ normal(beta[1],sigma);
}
'''
# Stan object 
StanDSO = ps.StanModel(model_code = model_string)

# data 
 data = dict(N = len(light_speed),
        y = light_speed
       )

#fit and verify model results 
fit_model = StanDSO.sampling(data=data, iter = 5000, chains = 2, thin = 1)
print fit_model

在pymc3

model_1= pm.Model()

with model_1:
    #priors as specified in stan model 
    mu = pm.Uniform('mu', lower = -np.inf, upper= np.inf)
    sigma = pm.Uniform('sigma', lower = 0, upper= np.inf)


    #using vague priors works 
    #mu = pm.Uniform('mu', lower = light_speed.std()/1000.0, upper= light_speed.std()*1000.0)
    #sigma = pm.Uniform('sigma', lower = light_speed.std()/1000.0, upper= light_speed.std()*1000.0)

    # define likelihood
    y_obs = pm.Normal('Y_obs', mu = mu, sd = sigma, observed = light_speed)

    # inference fitting the model

    # I have to use slice because the following command 
    #trace = pm.sample(5000)
    # produce the error 
    # ValueError: Cannot construct a ufunc with more than 32 operands 
    #(requested number were: inputs = 51 and outputs = 1)valueerror 


    xstart = pm.find_MAP()
    xstep = pm.Slice()
    trace = pm.sample(5000, xstep, xstart, random_seed = 123, progressbar= True)

pm.summary(trace)

我查看了 glm 的源代码，它似乎默认使用模糊的先验，Krutschke 和许多 BUGS 示例也推荐这种做法并且它有效（见上文）。然而，立场参考手册（第 52 页）表示，他们积极劝阻用户不要使用默认比例先验，因为他们将太多的概率质量集中在合理的后验值之外（...）可能会产生倾斜后验的深远影响。

4

1 回答 1

3

统一先验在 Stan 中定义在参数的支持上。因此，如果您声明一个参数，该参数real<lower=0> sigma;声明sigma对正值具有统一的先验。它通过对数转换sigma为 (-inf, inf) 然后考虑雅可比来做到这一点；我相信 PyMC3 可以做同样的事情。如果后验是正确的，Stan 允许不正确的先验。除了为变量指定的默认统一先验之外的任何其他先验都会相乘（在对数尺度上添加），因此它们的行为与预期一样（默认统一分布无效）。

话虽如此，我们建议在 Stan 手册中使用至少信息量较弱的先验，如果不是更强的先验。查看 Gelman 在回归部分的先验讨论中引用的论文，看看先验是如何被扭曲的。不仅仅是先验过于模糊；通常封闭区间上的统一先验（如在许多 BUGS 示例中）会对后验产生强烈影响——您可以根据截断似然函数来形象化这一点。所有这些都可能导致计算问题，具体取决于模型的后验几何形状（不仅适用于 HMC，也适用于 Gibbs 或 Metropolis）。

我们在 GitHub 上的 Gelman 和 Hill 回归书示例 (ARM) 不符合 Stan 当前的建模标准——它们只是直接从 ARM 代码翻译而来。修改所有这些都在我们 2016 年的待办事项清单上（以及 BUGS 示例）。一些 ARM 示例是从书中的 Rlm()或glm()函数示例翻译而来的，还有一些是lme4包lmer()和glmer()函数的示例，其中没有一个接受先验。我们即将发布 RStanARM 包，它将接受 R 的回归（和 lme4 的多级回归）表示法，并允许 MLE（有错误）或 HMC 估计（在后端具有高度优化的模型代码）。

于 2016-01-01T20:05:55.220 回答