将WINBUGS模型转换为PyMC3

问题描述

我目前正在上贝叶斯统计课程，我们可以使用任何软件包以计算方式求解模型，但是WINBUGS中提供了所有示例。我更喜欢使用Python和PyMC3。我在使用PyMC3方面没有太多经验，可以在如何将这种简单的WINBUGS模型转换为PyMC3模型方面提供一些帮助。

下面是示例WINBUGS代码。这是一个简单的二项式，将两个选项与每个样本的观察值数量进行比较。该模型还测试了5种不同的先验条件。

model{
for(i in 1:5){
   n1[i] <- Tot1   #100
   n2[i] <- Tot2   #3
   y1[i] <- Positives1
   y2[i] <- Positives2
   y1[i] ~ dbin(p1[i],n1[i])
   y2[i] ~ dbin(p2[i],n2[i])
  diffps[i] <- p1[i]-p2[i]   #100seller - 3seller
}

# Uniform priors
   p1[1] ~ dbeta(1,1); p2[1] ~ dbeta(1,1)

# Jeffreys'  priors
   p1[2] ~ dbeta(0.5,0.5);  p2[2] ~ dbeta(0.5,0.5)

# informative priors centered at about 93% and 97%
   p1[3] ~ dbeta(30,2);  p2[3] ~ dbeta(2.9,0.1)

# Zellner priors prop to 1/(p * (1-p))
  logit(p1[4]) <- x[1]
    x[1] ~ dunif(-10000,10000)  # as  dflat()
   logit(p2[4]) <- x[2]
   x[2] ~ dunif(-10000,10000)   # as dflat()

#Logit centered at 3 gives mean probs close to 95%
  logit(p1[5]) <- x[3]
    x[3] ~  dnorm(3,1)  
   logit(p2[5]) <- x[4]
    x[4] ~ dnorm(3,1) 
}

DATA
list(Tot1=100,Tot2=3,Positives1=95,Positives2=3)

INITS
list(p1=c(0.9,0.9,NA,NA),p2=c(0.9,x=c(0,0))

list(p1=c(0.5,0.5,p2=c(0.5,0))

list(p1=c(0.3,0.3,p2=c(0.3,0))

在PyMC3中，我尝试使用以下代码在单个样本上实现5个先验中的第一个（我不确定如何同时实现）：

import np as np
import pymc3 as pm

sample2 = np.ones(3)

with pm.Model() as ebay_example:
    prior = pm.Beta('theta',alpha = 1,beta = 1)
    likelihood = pm.Bernoulli('y',p = prior,observed = sample2)
    trace = pm.sample(1000,tune = 2000,target_accept = 0.95)

运行了上述模型，但是结果与BUGS结果不一致，我不确定是否是因为我没有进行老化或其他较大的问题。任何指导都是很棒的。

解决方法

我们现在正在上相同的课。以下是我的代码，均值差异接近BUGS结果。

na = 100
nb = 3
pos_a = 95
pos_b = 3

with pm.Model() as model:
    # priors
    p0a = pm.Beta('p0a',1,1)

    # likelihood
    obs_a = pm.Binomial("obs_a",n=na,p=p0a,observed=pos_a)

    # sample
    trace1_a = pm.sample(1000)


with pm.Model() as model:
    # priors
    p0b = pm.Beta('p0b',1)

    # likelihood
    obs_b = pm.Binomial("obs_b",n=nb,p=p0b,observed=pos_b)

    # sample
    trace1_b = pm.sample(1000)

pm.summary(trace1_a)["mean"][0] - pm.summary(trace1_b)["mean"][0]
OUT: 0.1409999999999999

pymc3 python winbugs