<jupyter_start><jupyter_text>Custom likelihoodsSometimes we want to be able to perform Bayesian inference on models which have not been implemented in probabilistic programming languages like `pymc3`. In these cases, we need to be able to sample from models for which we can write down the analytic log-likelihood ourselves.Tip: See `distribution.py` and `dist_math.py` in `~/anaconda3/lib/python3.6/site-packages/pymc3/distributions` for examples of how to implement custom likelihoods.Author: <NAME><jupyter_code>import pymc3 as pm
import numpy as np
import matplotlib.pyplot as plt
import mystyle.sty as sty # see https://github.com/jaryaman/mystyle 
sty.reset_plots()
from pymc3.distributions.dist_math import binomln, logpow
import theano.tensor as tt<jupyter_output><empty_output><jupyter_text>Suppose that we would like to sample from the following likelihood, which is the solution for a pure-birth process:$$P(i|\lambda t, n_0) = \begin{pmatrix}i-1\\n_0 - 1\end{pmatrix} e^{-\lambda n_0 t} (1- e^{-\lambda t})^{i - n_0}$$which is a negative binomial distribution. Define the log-likelihood for the model:<jupyter_code>def logp_neg_binom(n0, lam_t, i):
    '''
    Log-likelihood for pure-birth process
    '''
    return binomln(i-1, n0-1) -lam_t*n0 + logpow(1 - tt.exp(-lam_t), i - n0)<jupyter_output><empty_output><jupyter_text>Given values for $n_0$ and $i$, compute the maximum likelihood estimate for $\lambda t$. Note that the maximum likelihood is equivalent to the maximum a posteriori estimate given a uniform prior for the parameter, and also assuming that the prior has support over the maximum likelihood estimate. Compile the model<jupyter_code>n0 = 50
i = 100
with pm.Model() as neg_binom_model:
    lam_t = pm.Uniform('lam-t', 0, 100)
    ll = pm.DensityDist('ll', logp_neg_binom, observed=dict(n0=n0,lam_t=lam_t,i=i))<jupyter_output><empty_output><jupyter_text>Find the MAP<jupyter_code>map_estimate = pm.find_MAP(model=neg_binom_model)
map_est = map_estimate['lam-t']
map_est
n0*np.exp(map_est.flatten()[0])<jupyter_output><empty_output><jupyter_text>So the maximum likelihood estimate recovers the true value of `i` quite accurately in this case We can also do full Bayesian inference:<jupyter_code>with neg_binom_model:
    trace = pm.sample(5000)
pm.plots.autocorrplot(trace);
fig, axs = plt.subplots(1,2, figsize=(2*5,5))
pm.plots.traceplot(trace, ax=np.reshape(axs, (1,2)));
axs[0].plot(map_est.flatten()[0]*np.ones(50), np.linspace(0,4.3),'--k', label='ML')
axs[0].legend();<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11632492113564669, 'fraction_numerical': 0.02168769716088328, 'mean_word_length': 5.801608579088472, 'pattern_counts': {'":': 0, '<': 22, '<?xml version=': 0, '>': 22, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '57570', 'license': 'MIT', 'n_tokens_mistral': 910, 'n_tokens_neox': 864, 'n_words': 251, 'path': 'Notebooks/custom_likelihoods.ipynb', 'repo_name': 'jaryaman/ML_demos', 'size': '2543'}