pymc3의 중첩 된로 짓 모델의 올바른 지정

Question

선택 항목 A, B 또는 C를 나타내는 3 가지 결과가 포함 된 중첩 된 로지스틱 회귀 모델에서 작업하고 있습니다. 첫 번째 수준은 A와 B 또는 C 중에서 선택하고 두 번째 수준은 선택 항목 레벨은 B와 C 사이의 선택을 나타냅니다. 일부 데이터를 작성한 코드는 아래에 있지만 모델을 올바르게 지정했는지 확실하지 않습니다. B 또는 C의 확률은 정의에 의해 B의 확률보다 높지만 아주 작은 표본 크기의 사후에서 샘플링 할 때 "BorC"는 B보다 작을 수 있습니다. 이러한 작은 표본 크기는 실제 데이터에 나타나지 않을 것입니다 나는 관심이있다. 그러나 이것이 일어난다는 사실은 나에게 내가 틀린 일을했다고 생각하게 만든다. 감사!

import numpy as np 
import pandas as pd 
import pymc3 as pm 
from scipy.special import logit 
import matplotlib.pyplot as plt 
from theano import shared 

import cPickle as pickle # python 2 

def hierarchical_normal(name, shape, mu=0.,cs=5.): 
    delta = pm.Normal('delta_{}'.format(name), 0., 1., shape=shape) 
    sigma = pm.HalfCauchy('sigma_{}'.format(name), cs) 

    return pm.Deterministic(name, mu + delta * sigma) 

NUTS_KWARGS = {'target_accept': 0.99} 
SEED = 4260026 # from random.org, for reproducibility 

np.random.seed(SEED) 

ndraws = 1000 

counts =[[19, 50, 37], 
     [21, 67, 55], 
     [11, 53, 38], 
     [17, 54, 45], 
     [24, 93, 66], 
     [27, 53, 70]] 

counts = pd.DataFrame(counts,columns=["A","B","C"]) 

counts["BorC"] = counts[["B","C"]].sum(axis=1) 
counts["n"] = counts[["A","B","C"]].sum(axis=1) 
print counts 

group = counts.index.values 
n_group = np.unique(group).size 

obs_B = counts.B.values 
obs_BorC = counts.BorC.values 
obs_n = counts.n.values 

obs_n_ = shared(obs_n) 

with pm.Model() as model: 
    beta0  = pm.Normal('beta0', 0.,sd=5.) 
    beta_group = hierarchical_normal('beta_group', n_group) 

    eta_BorC = beta0 + beta_group[group] 
    p_BorC  = pm.math.sigmoid(eta_BorC) 
    like_BorC = pm.Binomial('obs_BorC', obs_n_, p_BorC, observed=obs_BorC) 

    alpha0  = pm.Normal('alpha0', 0.,sd=5.) 
    alpha_group = hierarchical_normal('alpha_group', n_group) 

    eta_BgivenBorC = alpha0 + alpha_group[group] 
    p_BgivenBorC = pm.math.sigmoid(eta_BgivenBorC) 
    like_BgivenBorC = pm.Binomial('obs_BgivenBorC', obs_BorC, p_BgivenBorC, observed=obs_B) 

    p_B = p_BgivenBorC*p_BorC 
    like_B = pm.Binomial('obs_B', obs_n_, p_B, observed=obs_B) 

    trace = pm.sample(draws=ndraws, random_seed=SEED,nuts_kwargs=NUTS_KWARGS) 

obs_n_.set_value(np.array([10]*6)) 

pp_trace = pm.sample_ppc(trace, model=model) 
C = pp_trace['obs_BorC']-pp_trace['obs_B'] 
print np.min(np.min(C)) 
B = pp_trace['obs_B'] 
C = np.sum(C,axis=1) 
B = np.sum(B,axis=1) 
diff = C-B 
plt.figure() 
plt.hist(diff,50) 
plt.show() 

편집 : 나는 로그 우도가 자동으로 like_B 내 사양은 중복을 의미 모든 변수를 합산하는 pymc3 코드를 검색에서 참조하십시오. 이 경우 후방 샘플링을 위해 obs_BorC를 올바르게 설정하는 방법을 알아야합니다.

Answer 1

시도한 해결책이 있습니다. 문제가 해결되지 않으면 해결 방법으로 사용할 수 있습니다. obs_BorC가 각 반복을 재설정하는 사용자 정의 사후 샘플러를 만들었습니다.

import numpy as np 
import pandas as pd 
import pymc3 as pm 
from scipy.special import logit 
import matplotlib.pyplot as plt 
from theano import shared 
from collections import defaultdict 

from scipy.stats import norm 

def hierarchical_normal(name, shape, mu=0.,cs=5.): 
    delta = pm.Normal('delta_{}'.format(name), 0., 1., shape=shape) 
    sigma = pm.HalfCauchy('sigma_{}'.format(name), cs) 

    return pm.Deterministic(name, mu + delta * sigma) 

NUTS_KWARGS = {'target_accept': 0.99} 
SEED = 4260026 # from random.org, for reproducibility 

np.random.seed(SEED) 

ndraws = 1000 

counts =[[19, 50, 37], 
     [21, 67, 55], 
     [11, 53, 38], 
     [17, 54, 45], 
     [24, 93, 66], 
     [27, 53, 70]] 

counts = pd.DataFrame(counts,columns=["A","B","C"]) 

counts["BorC"] = counts[["B","C"]].sum(axis=1) 
counts["n"] = counts[["A","B","C"]].sum(axis=1) 
print counts 

group = counts.index.values 
n_group = np.unique(group).size 

obs_B = counts.B.values 
obs_BorC = counts.BorC.values 
obs_n = counts.n.values 

obs_n_ = shared(obs_n) 
obs_BorC_ = shared(obs_BorC) 

with pm.Model() as model: 
    beta0  = pm.Normal('beta0', 0.,sd=5.) 
    beta_group = hierarchical_normal('beta_group', n_group) 

    eta_BorC = beta0 + beta_group[group] 
    p_BorC  = pm.math.sigmoid(eta_BorC) 
    like_BorC = pm.Binomial('obs_BorC', obs_n_, p_BorC, observed=obs_BorC) 

    alpha0  = pm.Normal('alpha0', 0.,sd=5.) 
    alpha_group = hierarchical_normal('alpha_group', n_group) 

    eta_BgivenBorC = alpha0 + alpha_group[group] 
    p_BgivenBorC = pm.math.sigmoid(eta_BgivenBorC) 
    like_BgivenBorC = pm.Binomial('obs_BgivenBorC', obs_BorC_, p_BgivenBorC, observed=obs_B) 

    trace = pm.sample(draws=ndraws, random_seed=SEED,nuts_kwargs=NUTS_KWARGS) 

#plt.figure() 
#axs = pm.forestplot(trace,varnames=['beta0','beta_group','alpha0','alpha_group']) 
#plt.savefig("Forest.png") 
#plt.close() 

obs_n_.set_value(np.array([10000]*6)) 
indices = np.random.randint(0, len(trace), 1000) 
ppc = defaultdict(list) 
for idx in indices: 
    param = trace[idx] 
    n_BorC = model["obs_BorC"].distribution.random(point=param,size=None) 
    obs_BorC_.set_value(np.array(n_BorC)) 

    n_B = model["obs_BgivenBorC"].distribution.random(point=param,size=None) 
    ppc["obs_BorC"].append(n_BorC) 
    ppc["obs_B"].append(n_B) 
pp_trace = {k: np.asarray(v) for k, v in ppc.items()} 
#pp_trace = pm.sample_ppc(trace, model=model,samples=20000) 
C = pp_trace['obs_BorC']-pp_trace['obs_B'] 
print np.min(np.min(C)) 
B = pp_trace['obs_B'] 
C = np.sum(C,axis=1) 
B = np.sum(B,axis=1) 
diff = (C-B)/60000. 
plt.figure() 
x = np.linspace(np.min(diff),np.max(diff),200) 
mean = np.mean(diff) 
std = np.std(diff) 
plt.hist(diff,50,normed=True) 
plt.plot(x,norm.pdf(x,mean,std)) 
plt.plot() 
plt.show()