# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python [conda env:root] *
#     language: python
#     name: conda-root-py
# ---

# + [markdown] slideshow={"slide_type": "slide"}
# ## PyMC3 example

# +
import numpy as np
import pandas as pd
import pymc3 as pm

# %config InlineBackend.figure_format = 'retina'
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()
# -

import sys
print("Python version:", sys.version)
print("Numpy version:", np.__version__)
print("PyMC3 version:", pm.__version__)

df = pd.read_csv("intervals.csv")

df.head()

# ![](reich_figure_1.png)
#
# Figure 1 from Reich et al. (2009), _Estimating incubation period distributions with coarse data_. 

Tmin = np.array(np.maximum(df["SL"]-df["ER"], 0))
Tmax = np.array(df["SR"]-df["EL"])

plt.hist(Tmin);

plt.hist(Tmax);

# +
xs = np.linspace(0, 75, 500)
ys = np.zeros(len(xs))

nCases = len(Tmin)
for i in range(nCases):
    ys[(xs >= Tmin[i]) & (xs <= Tmax[i])] += 1 / (nCases * (Tmax[i]-Tmin[i]))
    
plt.plot(xs, ys);
# -

Tavgs = (Tmin + Tmax) / 2

# %%time
with pm.Model() as model:
    μ = pm.Uniform('μ', lower=-25, upper=25)
    σ = pm.Uniform('σ', lower=0, upper=25) 
    T = pm.Lognormal('T', mu=μ, sigma=σ, observed=Tavgs)
    trace = pm.sample(10**4)

pm.plot_trace(trace);

pm.plot_posterior(trace);

# This is assuming we have observations $\boldsymbol{T} = (T_1, \dots, T_n)$ which gives us a likelihood of
#
# $$ L(\mu, \sigma \mid \boldsymbol{T}) = \prod_{i=1}^n f(T_i; \mu, \sigma )  $$
#
# where $f(x; \mu, \sigma)$ is the p.d.f. of the $\mathsf{LogNormal}(\mu, \sigma^2)$ distribution.
#
# However we don't have observations, just intervals. Say each unobserved period fell into $T_i \in [T_i^-, T_i^+]$. Our likelihood becomes
#
# $$ L(\mu, \sigma \mid \boldsymbol{T}^{\,-}, \boldsymbol{T}^{\,+}) = \prod_{i=1}^n \bigl[ F(T_i^+ ; \mu, \sigma) - F(T_i^-; \mu, \sigma) \bigr] $$
#
# where $F(x; \mu, \sigma)$ is the c.d.f. of the $\mathsf{LogNormal}(\mu, \sigma^2)$ distribution.

# +
import theano.tensor as tt

# Taken from PyMC3's pymc3/distributions/dist_math.py file
# starting at line 346.
def zvalue(x, sigma, mu):
    """
    Calculate the z-value for a normal distribution.
    """
    return (x - mu) / sigma

# Taken from PyMC3's pymc3/distributions/continuous.py file
# starting at line 1849.
def cdf(x, mu, sigma):
        """
        Compute the log of the cumulative distribution function for Lognormal distribution
        at the specified value.

        Parameters
        ----------
        x: numeric
            Value(s) for which log CDF is calculated. If the log CDF for multiple
            values are desired the values must be provided in a numpy array or theano tensor.

        Returns
        -------
        TensorVariable
        """
        z = zvalue(np.log(x), mu=mu, sigma=sigma)

        return tt.switch(
                tt.lt(z, -1.0),
                tt.erfcx(-z / tt.sqrt(2.)) / 2. * np.exp(-tt.sqr(z) / 2),
                tt.erfc(-z / tt.sqrt(2.)) / 2.
        )


# -

# With `Potential` we have to add log-terms to the likelihood. So
#
#
# $$ \log \bigl[ L(\mu, \sigma \mid \boldsymbol{T}^{\,-}, \boldsymbol{T}^{\,+}) \bigr] = \sum_{i=1}^n \log \bigl[  F(T_i^+ ; \mu, \sigma) - F(T_i^-; \mu, \sigma) \bigr] \,. $$
#

# %%time
with pm.Model() as model:
    μ = pm.Uniform('μ', lower=-25, upper=25)
    σ = pm.Uniform('σ', lower=0, upper=25) 
    pm.Potential('T', tt.sum(tt.log( cdf(Tmax, μ, σ) - cdf(Tmin, μ, σ) )))
    trace = pm.sample(10**5, step=pm.Metropolis()) 

pm.stats.ess(trace["μ"]), pm.stats.ess(trace["σ"]) 

pm.plot_trace(trace);

pm.plot_posterior(trace);

trace["μ"].mean()

trace["σ"].mean()

StarcoderdataJupyterScriptsDedupFiltered