I am new to PyMC1, and thought I understood it 
, but obviously not.
In the below, I encode a simple discrete Bayesian model with no hyper-parameters. I encode the observations in an observation, which is passed to the model. In the experiments below, I set ‘l’ to 0 and to ‘1’. ‘l’ is one of the two possible results from sampling a Bernoulli. The ‘l’:0 case does not have a problem. The ‘l’:1 case does. I read up on this error, and find that it often occurs when a parameter is out of range with respects to constrains imposed by some distribution (for example, a negative standard deviation). However, in the present case, I cannot figure it out. Any help is very much appreciated. Thanks. Gordon.
import pymc3 as pm
import pymc3.distributions.discrete as discrete
import matplotlib.pyplot as plt
import arviz as az
# import aesara.tensor as at
import numpy as np
from tqdm import tqdm
import pandas as pd
from theano import tensor as tt
from theano.ifelse import ifelse
# SOME DETERMINISTIC FUNCTIONS
def grade(i, d):
a00 = tt.and_(tt.eq(i,0), tt.eq(d,0))
a10 = tt.and_(tt.eq(i,1), tt.eq(d,0))
a01 = tt.and_(tt.eq(i,0), tt.eq(d,1))
a11 = tt.and_(tt.eq(i,1), tt.eq(d,1))
cat1 = [0.3,0.4, 0.3]
cat2 = [0.05, 0.25, 0.7]
cat3 = [0.9, 0.08, 0.02]
cat4 = [0.5, 0.3, 0.2]
return a00*cat1 + a10*cat2 + a01*cat3 + a11*cat4
def letter(g):
return tt.eq(g, gA)*0.9 + tt.eq(g, gB)*0.6 + tt.eq(g, gC)*0.01
def score(i):
return tt.eq(i, 0)*0.05 + tt.eq(i, 1)*0.8
gA = 2
gB = 5
gC = 10
# MODEL CREATION: observations passed in as arguments through a dictionary
def create_model(obs):
with pm.Model() as grade_model:
if "d" in obs.keys():
d = discrete.Bernoulli("d", 0.4, observed=obs["d"])
else:
d = discrete.Bernoulli("d", 0.4)
if "i" in obs.keys():
i = discrete.Bernoulli("i", 0.3, observed=obs["i"])
else:
i = discrete.Bernoulli("i", 0.3)
p_grade = grade(i, d)
if "g" in obs.keys():
g = discrete.Categorical("g", p_grade, observed=obs["g"])
else:
g = discrete.Categorical("g", p_grade)
p_letter = letter(g)
if "l" in obs.keys():
l = discrete.Bernoulli("l", p_letter, observed=obs["l"])
else:
l = discrete.Bernoulli("l", p_letter)
print("l: ", l)
p_score = score(i)
if "s" in obs.keys():
s = discrete.Bernoulli("s", p_score, observed=obs["s"])
else:
s = discrete.Bernoulli("s", p_score)
return grade_model;
#
def run_model(model):
with model:
step = pm.Metropolis()
trace = pm.sample(n_init=5000, step=step, draws=5000, tune=1000, return_inferencedata=False, progressbar=False)
display(az.summary(trace, kind="stats"))
lst = list(trace.points())
# glst = [dct['g'] for dct in lst]
# slst = [dct['s'] for dct in lst]
# llst = [dct['l'] for dct in lst]
return lst, trace #, [np.mean(x) for x in [glst, slst, llst]];
def run(obs):
model = create_model(obs)
# model = create_model_new(obs)
lst, trace = run_model(model)
return model, trace
def run(obs):
model = create_model(obs)
# model = create_model_new(obs)
lst, trace = run_model(model)
return model, trace
model1, trace1 = run(obs={})
# RUNS FINE.
run(obs={'l':0.})
# PRODUCES ERROR (the observation is in a Bernoulli)
run(obs='l':1}) # <<< ERROR
This last line produces the following error:
---------------------------------------------------------------------------
SamplingError Traceback (most recent call last)
<ipython-input-12-93d95129a24b> in <module>
1 # 'l':1 gives an error related to starting conditions
----> 2 run(obs={'l':1})
<ipython-input-8-364b5c5a1e7d> in run(obs)
2 model = create_model(obs)
3 # model = create_model_new(obs)
----> 4 lst, trace = run_model(model)
5 return model, trace
<ipython-input-7-f445ae35680c> in run_model(model)
3 # return_inferencedata=True does not work
4 step = pm.Metropolis()
----> 5 trace = pm.sample(n_init=5000, step=step, draws=5000, tune=1000, return_inferencedata=False, progressbar=False)
6 display(az.summary(trace, kind="stats"))
7 lst = list(trace.points())
~/anaconda3/lib/python3.8/site-packages/deprecat/classic.py in wrapper_function(wrapped_, instance_, args_, kwargs_)
213 else:
214 warnings.warn(message, category=category, stacklevel=_routine_stacklevel)
--> 215 return wrapped_(*args_, **kwargs_)
216
217 return wrapper_function(wrapped)
~/anaconda3/lib/python3.8/site-packages/pymc3/sampling.py in sample(draws, step, init, n_init, initvals, trace, chain_idx, chains, cores, tune, progressbar, model, random_seed, discard_tuned_samples, compute_convergence_checks, callback, jitter_max_retries, start, return_inferencedata, idata_kwargs, mp_ctx, pickle_backend, **kwargs)
442 start = deepcopy(start)
443 if start is None:
--> 444 check_start_vals(model.test_point, model)
445 else:
446 if isinstance(start, dict):
~/anaconda3/lib/python3.8/site-packages/pymc3/util.py in check_start_vals(start, model)
235
236 if not np.all(np.isfinite(initial_eval)):
--> 237 raise SamplingError(
238 "Initial evaluation of model at starting point failed!\n"
239 "Starting values:\n{}\n\n"
SamplingError: Initial evaluation of model at starting point failed!
Starting values:
{'d': array(0), 'i': array(0), 'g': array(1), 's': array(0)}
Initial evaluation results:
d -0.51
i -0.36
g -0.92
s -0.05
l -inf
Name: Log-probability of test_point, dtype: float64