import pymc3 as pm
import matplotlib.pyplot as plt
import numpy as np
from utils import ECDF
from data import load_decay
import pandas as pd
import theano.tensor as tt
import arviz as az
%load_ext autoreload
%autoreload 2
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
Now that you've learned about Bayesian estimation, we're going to explore one more topic: Bayesian curve fitting.
By "curve fitting", we're really talking about any curve: those that are bendy, those that are straight, and those that are in between.
In order to reinforce this point, rather than show you plain vanilla linear regression, we will work through an exponential decay curve example.
You've taken radioactive decay measurements of an unknown element in a secure facility. The measurements are noisy, though, and potentially have some bias. In the face of this, we would like to be able to characterize the decay constant of this unknown material, potentially leading to an identification of the material.
Let's load in the data.
np.random.seed(42)
df = load_decay()
df.head(5)
t | activity | |
---|---|---|
0 | 0 | 63.496714 |
1 | 1 | 62.281634 |
2 | 2 | 62.495498 |
3 | 3 | 62.806653 |
4 | 4 | 60.493075 |
Plot activity
vs. time
.
ax = df['activity'].plot()
If we were to draw out a model for the curve above, how might it look like?
(To reveal one possible model, double-click on this Markdown cell and remove the z
from the end of the filename.)
The most important part of this diagram is the "link function" - this is what "links" the data to the output. In this case, we've used the exponential decay curve as the link function, but if you were doing a linear regression model, all you would have to do is to change the link function for the $y=mx+c$ "straight curve", and do another curve fit with the appropriate priors for $m$ and $c$.
If you're familiar with the mathematical groundings of deep learning, you'll immediately recognize that a deep neural network model is merely another instance of a really complicated link function that links the input data $x$ to the observed data $y$, with the model weights and biases corresponding to the parameters (let's collectively call this set of parameters $\theta$).
Now that you've seen a pictorial description of the model, implement it below in PyMC3.
with pm.Model() as model:
A = pm.HalfNormal('A', sd=100)
tau = pm.Exponential('tau', lam=1)
C = pm.Normal('C', sd=100)
sd = pm.HalfCauchy('sd', beta=1)
link = A * np.exp(-df['t'].values / tau) + C
like = pm.Normal('activity', mu=link, sd=sd, observed=df['activity'].values)
with model:
trace = pm.sample(2000, tune=2000)
# Note: Sampler may pause for a while after finishing
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (4 chains in 4 jobs) NUTS: [sd, C, tau, A] Sampling 4 chains: 100%|██████████| 16000/16000 [00:07<00:00, 2272.16draws/s]
Check that sampling has converged.
az.plot_trace(trace);
az.plot_posterior(trace);
More generally, if
$$y = f(x, \theta)$$where $\theta$ are merely a set of parameters, then you can perform inference on the curve's parameters $\theta$. To make this clear:
curve name | functional form | parameters |
---|---|---|
exponential decay | $y = Ae^{-t/\tau} + C$ | $A$, $\tau$, $C$ |
sine curves | $y = A\sin(\omega x + \phi)$ | $A$, $\omega$, $\phi$ |
linear regression | $y = mx + c$ | $m$, $c$ |
logistic regression | $y = L(mx + c)$ | $m$, $c$ |
4-parameter IC50 | $y = \frac{a - i}{1 + 10^{\beta(log(\tau) - x)}} + i$ | $a$, $i$, $\tau$, $\beta$ |
deep learning | $y = f(x, \theta)$ | $\theta$ |