Indexing using a free random variable

Questions
15

I had another look at the model, with the following code:

with pm.Model() as basic_model:
    # Priors for unknown model parameters
    b1 = pm.Uniform('b1', lower=0.3, upper=0.5, testval=0.45)
    ncomp_aCH = pm.Categorical('ncomp_aCH', p=np.ones(n)/n)
    ncomp_aCOH = pm.Categorical('ncomp_aCOH', p=np.ones(n)/n)

    aCH=aCH_[0, ncomp_aCH]
    aCOH=aCOH_[0, ncomp_aCOH]

    out= b1*aCH+aCOH

    # Likelihood (sampling distribution) of observations
    Y_obs = pm.Normal('Y_obs', mu=out, tau=sigma, observed=Y)
    trace = pm.sample(discard_tuned_samples=False)

which gives me a trace:

image
image.png856×424 37.9 KB

So looking at the trace, The problem seems to be that ncomp_aCH is stuck after a few step. Either it is highly likely (think of the posterior being almost like a delta function) or there is some problem with the sampler.

I verify the logp of ncomp_aCH using the following:

point = trace[-1]
point
# {'b1': 0.4980582094475776,
#  'b1_interval__': 4.624950462026797,
#  'ncomp_aCH': 14,
#  'ncomp_aCOH': 17}
for i in range(n):
    point['ncomp_aCH']=i
    print(point['ncomp_aCH'], basic_model.logp(point))

0 -1110.9799301118057
1 -1118.7101039762993
2 -1128.3039141889512
3 -1125.591418583891
4 -1133.5723103615805
5 -1134.252739856408
6 -1131.6183448546024
7 -1132.7872247881817
8 -1111.4298570895123
9 -1134.4647708083562
10 -1144.1176821443162
11 -1141.1537636543512
12 -1097.8710713143184
13 -1086.384805832985
14 -1038.9504119754274
15 -1053.7949113856769
16 -1102.021203844887
17 -1102.374012701112
18 -1104.585462412219
19 -1148.0232245001032

Looking from the output, ncomp_aCH=14 does have the lowest logp. In fact, the dlogp between ncomp_aCH=14 and ncomp_aCH=15 is np.exp(1038.9504119754274-1053.7949113856769) = 3.573681299085053e-07 (the acceptance proability), which is why all the proposal value are rejected.
I think you should double check your model, as the output from the sampler seems to be correct.