^ I will be discussing general implementations of Bayesian methods
... the Bayesian approach is attractive because it is useful. Its usefulness derives in large measure from its simplicity. Its simplicity allows the investigation of far more complex models than can be handled by the tools in the classical toolbox.
Link and Barker 2010
^ worth the trouble!
^ Sampling from un-normalized distributions
Markov chain Monte Carlo simulates a
for which some function of interest is the
distribution.
^ e.g. the joint distribution of the parameters of some model
This is guaranteed when the Markov chain is constructed that satisfies the detailed balance equation:
^ First-generation MCMC
model {
for (j in 1:J){
y[j] ~ dnorm (theta[j], tau.y[j])
theta[j] ~ dnorm (mu.theta, tau.theta)
tau.y[j] <- pow(sigma.y[j], -2)
}
mu.theta ~ dnorm (0.0, 1.0E-6)
tau.theta <- pow(sigma.theta, -2)
sigma.theta ~ dunif (0, 1000)
}
^
^ Debugging is impossible
^ Andrew Thomas
OpenBUGS is open source only in the read-only sense. -- Andrew Thomas
^ Source code is not posted on public website
Atrocities include:
^
A platform for experimenting with ideas in Bayesian modeling
^
proc mcmc data=stagnant outpost=postout seed=24860 ntu=1000
nmc=20000;
ods select PostSummaries;
ods output PostSummaries=ds;
array beta[2];
parms alpha cp beta1 beta2;
parms s2;
prior cp ~ unif(-1.3, 1.1);
prior s2 ~ uniform(0, 5);
prior alpha beta: ~ normal(0, v = 1e6);
j = 1 + (x >= cp);
mu = alpha + beta[j] * (x - cp);
model y ~ normal(mu, var=s2);
run;
^
Gelman Lab, Columbia University
^ Stanislaw Ulam (& Nicholas Metropolis) -- Manhattan project
data {
int<lower=0> J; // number of schools
real y[J]; // estimated treatment effects
real<lower=0> sigma[J]; // s.e. of effect estimates
}
parameters {
real mu;
real<lower=0> tau;
real eta[J];
}
^
transformed parameters {
real theta[J];
for (j in 1:J)
theta[j] <- mu + tau * eta[j];
}
model {
eta ~ normal(0, 1);
y ~ normal(theta, sigma);
}
>>> theta = Exponential('theta', 2.)
>>> theta
<pymc.distributions.Exponential 'theta' at 0x10f9e2790>
>>>
>>> theta = Exponential('theta', 2.)
>>> theta
<pymc.distributions.Exponential 'theta' at 0x10f9e2790>
>>> theta.value
array(0.3848267102336743)
>>>
>>> theta = Exponential('theta', 2.)
>>> theta
<pymc.distributions.Exponential 'theta' at 0x10f9e2790>
>>> theta.value
array(0.3848267102336743)
>>> theta.logp
-0.07650623990740335
>>>
>>> theta = Exponential('theta', 2.)
>>> theta
<pymc.distributions.Exponential 'theta' at 0x10f9e2790>
>>> theta.value
array(0.3848267102336743)
>>> theta.logp
-0.07650623990740335
>>> theta.random()
array(0.4290156681214033)
>>>
>>> theta = Exponential('theta', 2.)
>>> theta
<pymc.distributions.Exponential 'theta' at 0x10f9e2790>
>>> theta.value
array(0.3848267102336743)
>>> theta.logp
-0.07650623990740335
>>> theta.random()
array(0.4290156681214033)
>>> theta.value
array(0.4290156681214033)
>>>
>>> theta = Exponential('theta', 2.)
>>> theta
<pymc.distributions.Exponential 'theta' at 0x10f9e2790>
>>> theta.value
array(0.3848267102336743)
>>> theta.logp
-0.07650623990740335
>>> theta.random()
array(0.4290156681214033)
>>> theta.value
array(0.4290156681214033)
>>> theta.logp
-0.1648841556828613
@observed
def survival(value=t, lam=theta, f=failure):
"""Exponential survival likelihood,
accounting for censoring"""
return sum(f * log(lam) - lam * value)
Christopher J. Fonnesbeck Melissa L. McPheeters Shanthi Krishnaswami Mary Louise Lindegren Tyler Reimschisel
J Inherit Metab Dis DOI 10.1007/s10545-012-9564-0
^part of a larger systematic review on the comparative effectiveness of treatment for PKU by Vanderbilt EPC
^ metabolic disorder in which a buildup of phenylalanine (Phe) in the blood results from an inability to properly metabolize it
^ Complete re-write of PyMC2
Amazon
Brown University
Quantopian
Uses a physical analogy of a frictionless particle moving on a hyper-surface
Requires an auxiliary variable to be specified
^
The Hamiltonian is defined as the sum of potential energy E(s) and kinetic energy K(φ)
^
^
with pm.Model() as model:
pm.glm.glm('y ~ x', data)
step = pm.Metropolis()
iter_sample = pm.iter_sample(samples+1000, step)
animation.FuncAnimation(fig, animate, init_func=init,
frames=samples, interval=5, blit=True)
with pm.Model() as model:
pm.glm.glm('y ~ x', data)
step = pm.NUTS()
iter_sample = pm.iter_sample(samples+1000, step)
animation.FuncAnimation(fig, animate, init_func=init,
frames=samples, interval=5, blit=True)
Hamiltonian MC requires tuning of two parameters:
^
NUTS uses a recursive algorithm to build a set of likely candidate points that spans a wide swath of the target distribution.
It stops automatically when it starts to double back and retrace its steps.
Doubling process builds a balanced binary tree whose leaf nodes correspond to position-momentum states.
Doubling is halted when the subtrajectory from the leftmost to the rightmost nodes of any balanced subtree of the overall binary tree starts to double back on itself
^
(Hoffmann and Gelman 2011)
^
200-dimensional mvn with high correlation
GOING TO HAVE A BAD TIME
>>> from theano import function,tensor as T
>>> x = T.dmatrix('x')
>>> s = T.sum(1 / (1 + T.exp(-x)))
>>> gs = T.grad(s, x)
>>> dlogistic = function([x], gs)
>>> dlogistic([[3, -1],[0, 2]])
array([[ 0.04517666, 0.19661193],
[ 0.25 , 0.10499359]])
^
Very little competition (or diversity) in space of Bayesian statistical software
Has potential for being "platform for experimenting with Bayesian ideas"
Much unfinished work