Using BetaBinomial in PyMC3
<p>I have a table of counts of binary outcomes and I would like to fit a beta binomial distribution to estimate $\alpha$ and $\beta$ parameters, but I am getting errors when I try to fit/sample the model distribution the way I do for other cases:</p>

<pre><code>import pymc3 as pm
import pandas as pd

df = pd.read_csv('~/data.csv', low_memory=False)
df = df[df.Clicks &gt;= 0]

C0=df.C.values
I0=df.N.values
N0 = C0 + I0

with pm.Model() as model:
    C=pm.constant(C0)
    I=pm.constant(I0)
    C1=pm.constant(C0 + 1)
    I1=pm.constant(I0 + 1)
    N=pm.constant(N0)
    alpha = pm.Exponential('alpha', 1/(C0.sum()+1))
    beta = pm.Exponential('beta', 1/(I0.sum()+1))
    obs = pm.BetaBinomial('obs', alpha, beta, N, observed=C0)


with model:
    advi_fit = pm.variational.advi(n=int(1e4))
    trace1 = pm.variational.sample_vp(advi_fit, draws=int(1e4))

pm.traceplot(trace1[::10])

with model:
    step = pm.NUTS()
    #step = pm.Metropolis() # &lt;== same problem
    trace2 = pm.sample(int(1e3), step)

pm.traceplot(trace2[::10])
</code></pre>

<p>In both cases the sampling fails with:</p>

<pre><code>MissingInputError: ("An input of the graph, used to compute Elemwise{neg,no_inplace}(P_logodds_), was not provided and not given a value.Use the Theano flag exception_verbosity='high',for more information on this error.", P_logodds
</code></pre>

<p>In the <code>advi</code> case the full stack trace is:</p>

<hr>

<pre><code>MissingInputError                         Traceback (most recent call last)
&lt;ipython-input-46-8947c7c798e5&gt; in &lt;module&gt;()
----&gt; 1 import codecs, os;__pyfile = codecs.open('''/tmp/py7996Jip''', encoding='''utf-8''');__code = __pyfile.read().encode('''utf-8''');__pyfile.close();os.remove('''/tmp/py7996Jip''');exec(compile(__code, '''/home/dmahler/Scripts/adops-bayes2.py''', 'exec'));

/home/dmahler/Scripts/adops-bayes2.py in &lt;module&gt;()
     59     advi_fit = pm.variational.advi(n=int(J*6.4e4), learning_rate=1e-3/J, epsilon=1e-8, accurate_elbo=False)
     60     #advi_fit = pm.variational.advi_minibatch(minibatch_RVs=[alpha, beta, p], minibatch_tensors=[C,I,N])
---&gt; 61     trace = pm.variational.sample_vp(advi_fit, draws=int(2e4))
     62 
     63 pm.traceplot(trace[::10])

/home/dmahler/Anaconda/lib/python2.7/site-packages/pymc3/variational/advi.pyc in sample_vp(vparams, draws, model, local_RVs, random_seed, hide_transformed)
    317 
    318     varnames = [str(var) for var in model.unobserved_RVs]
--&gt; 319     trace = NDArray(model=model, vars=vars_sampled)
    320     trace.setup(draws=draws, chain=0)
    321 

/home/dmahler/Anaconda/lib/python2.7/site-packages/pymc3/backends/ndarray.pyc in __init__(self, name, model, vars)
     21     """
     22     def __init__(self, name=None, model=None, vars=None):
---&gt; 23         super(NDArray, self).__init__(name, model, vars)
     24         self.draw_idx = 0
     25         self.draws = None

/home/dmahler/Anaconda/lib/python2.7/site-packages/pymc3/backends/base.pyc in __init__(self, name, model, vars)
     34         self.vars = vars
     35         self.varnames = [var.name for var in vars]
---&gt; 36         self.fn = model.fastfn(vars)
     37 
     38 

/home/dmahler/Anaconda/lib/python2.7/site-packages/pymc3/model.pyc in fastfn(self, outs, mode, *args, **kwargs)
    374         Compiled Theano function as point function.
    375         """
--&gt; 376         f = self.makefn(outs, mode, *args, **kwargs)
    377         return FastPointFunc(f)
    378 

/home/dmahler/Anaconda/lib/python2.7/site-packages/pymc3/memoize.pyc in memoizer(*args, **kwargs)
     12 
     13         if key not in cache:
---&gt; 14             cache[key] = obj(*args, **kwargs)
     15 
     16         return cache[key]

/home/dmahler/Anaconda/lib/python2.7/site-packages/pymc3/model.pyc in makefn(self, outs, mode, *args, **kwargs)
    344                                on_unused_input='ignore',
    345                                accept_inplace=True,
--&gt; 346                                mode=mode, *args, **kwargs)
    347 
    348     def fn(self, outs, mode=None, *args, **kwargs):

/home/dmahler/Anaconda/lib/python2.7/site-packages/theano/compile/function.pyc in function(inputs, outputs, mode, updates, givens, no_default_updates, accept_inplace, name, rebuild_strict, allow_input_downcast, profile, on_unused_input)
    318                    on_unused_input=on_unused_input,
    319                    profile=profile,
--&gt; 320                    output_keys=output_keys)
    321     # We need to add the flag check_aliased inputs if we have any mutable or
    322     # borrowed used defined inputs

/home/dmahler/Anaconda/lib/python2.7/site-packages/theano/compile/pfunc.pyc in pfunc(params, outputs, mode, updates, givens, no_default_updates, accept_inplace, name, rebuild_strict, allow_input_downcast, profile, on_unused_input, output_keys)
    477                          accept_inplace=accept_inplace, name=name,
    478                          profile=profile, on_unused_input=on_unused_input,
--&gt; 479                          output_keys=output_keys)
    480 
    481 

/home/dmahler/Anaconda/lib/python2.7/site-packages/theano/compile/function_module.pyc in orig_function(inputs, outputs, mode, accept_inplace, name, profile, on_unused_input, output_keys)
   1774                    profile=profile,
   1775                    on_unused_input=on_unused_input,
-&gt; 1776                    output_keys=output_keys).create(
   1777             defaults)
   1778 

/home/dmahler/Anaconda/lib/python2.7/site-packages/theano/compile/function_module.pyc in __init__(self, inputs, outputs, mode, accept_inplace, function_builder, profile, on_unused_input, fgraph, output_keys)
   1426             # OUTPUT VARIABLES)
   1427             fgraph, additional_outputs = std_fgraph(inputs, outputs,
-&gt; 1428                                                     accept_inplace)
   1429             fgraph.profile = profile
   1430         else:

/home/dmahler/Anaconda/lib/python2.7/site-packages/theano/compile/function_module.pyc in std_fgraph(input_specs, output_specs, accept_inplace)
    175 
    176     fgraph = gof.fg.FunctionGraph(orig_inputs, orig_outputs,
--&gt; 177                                   update_mapping=update_mapping)
    178 
    179     for node in fgraph.apply_nodes:

/home/dmahler/Anaconda/lib/python2.7/site-packages/theano/gof/fg.pyc in __init__(self, inputs, outputs, features, clone, update_mapping)
    169 
    170         for output in outputs:
--&gt; 171             self.__import_r__(output, reason="init")
    172         for i, output in enumerate(outputs):
    173             output.clients.append(('output', i))

/home/dmahler/Anaconda/lib/python2.7/site-packages/theano/gof/fg.pyc in __import_r__(self, variable, reason)
    358         # Imports the owners of the variables
    359         if variable.owner and variable.owner not in self.apply_nodes:
--&gt; 360                 self.__import__(variable.owner, reason=reason)
    361         if (variable.owner is None and
    362                 not isinstance(variable, graph.Constant) and

/home/dmahler/Anaconda/lib/python2.7/site-packages/theano/gof/fg.pyc in __import__(self, apply_node, check, reason)
    472                             "for more information on this error."
    473                             % str(node)),
--&gt; 474                             r)
    475 
    476         for node in new_nodes:

MissingInputError: ("An input of the graph, used to compute Elemwise{neg,no_inplace}(P_logodds_), was not provided and not given a value.Use the Theano flag exception_verbosity='high',for more information on this error.", P_logodds_)
&gt; /home/dmahler/Anaconda/lib/python2.7/site-packages/theano/gof/fg.py(474)__import__()
    472                             "for more information on this error."
    473                             % str(node)),
--&gt; 474                             r)
    475 
    476         for node in new_nodes:
</code></pre>

<p>Before I was aware of <code>pymc3.BetaBinomial</code> I was fitting trying to achieve the same result using separate <code>Beta</code> and <code>Binomial</code> instances:</p>

<pre><code>with pm.Model() as model:
    C=pm.constant(C0)
    I=pm.constant(I0)
    C1=pm.constant(C0 + 1)
    I1=pm.constant(I0 + 1)
    N=pm.constant(N0)
    alpha = pm.Exponential('alpha', 1/(C0.sum()+1))
    beta = pm.Exponential('beta', 1/(I0.sum()+1))
    p = pm.Beta('P', alpha, beta,  shape=K)
    b = pm.Binomial('B', N, p, observed=C0)
</code></pre>

<p>This completes successfully but different methods produce rather different results. I thought this could be partly due to the extra level of indirection between the priors and the observations makes the search space larger. When I came across <code>BetaBinomial</code> I figured it would make the search easier as well as being <em>the right thing ^TM</em>. Otherwise I believe the to models should be logically equivalent. Unfortunatelly I cannot figure out how to make <code>batebinomial</code> work and I am unable to find any examples using <code>BetaBinomial</code> on the internet.</p>

<ul>
<li>How do I make the <code>BetaBinomial</code> modwel work?</li>
<li>Are the models really logically equivalent?</li>
<li>Does anybody have a better guess at the cause of the numerical problems with the initial hierarchical version?

<ul>
<li>How I could fix them?</li>
</ul></li>
</ul>

<p>Your model should run, and you can write this </p>

<pre><code>with pm.Model() as model:
    alpha = pm.Exponential('alpha', 1/(C0.sum()+1))
    beta = pm.Exponential('beta', 1/(I0.sum()+1))
    obs = pm.BetaBinomial('obs', alpha, beta, N0, observed=C0)
</code></pre>

<p>That is, (C, I C1, I1) were defined in you model but not used. Anyway that was not the problem. The error you are getting is because PyMC3 is expecting a variable <code>P</code> (like in the second model you have) but the variable is not defined. May be you are working with a Jupyter notebook and you delete/comment a theano variable. Try running the Notebook again.</p>

<p>Theoretically using a Beta and a Binomial or a BetaBinomial should give the same results. From a practical point of view. Sampling from the BetaBinomial should be faster than from a beta and a binomial, since part of the work has already being done analytically!</p>

<p>Assuming a correct sampling both models should provide equivalent results. To check both results are equivalent try increasing the number of samples (and avoid <a href="http://doingbayesiandataanalysis.blogspot.co.at/2011/11/thinning-to-reduce-autocorrelation.html" rel="nofollow">thinning</a>). Also compare the results between and within models (the differences should be about the same). If you don't need to estimate the <code>P</code> variable (the beta distribution), then use the BetaBinomial.</p>

python bayesian pymc data-science pymc3