emcee.Originals available at: https://github.com/dfm/gp-tutorial & https://speakerdeck.com/dfm/an-astronomers-introduction-to-gaussian-processes
In [16]:
%pylab inline
# Auto imports np and plt
In [17]:
from matplotlib import rcParams
#rcParams["savefig.dpi"] = 150
import emcee # http://dan.iel.fm/emcee
import triangle # https://github.com/dfm/triangle.py
import numpy as np
import matplotlib.pyplot as plt
from astroML.plotting import hist
from astroML.stats.random import bivariate_normal
from matplotlib.patches import Ellipse
import timeit
#from IPython.display import display, Math, Latex
from astroML.decorators import pickle_results
In [18]:
import seaborn as sns
sns.set_context("notebook", font_scale=1.5, rc={"lines.linewidth": 2.5})
cmap = sns.cubehelix_palette(light=1, as_cmap=True)
sns.palplot(sns.cubehelix_palette(light=1))
In [19]:
import os
In [20]:
import pandas as pd
AO_dir = os.getenv("AO_bonding_paper")
df = pd.read_csv(AO_dir+"/data/cln_20130218_cary5000.csv")
df.set_index('wavelength', inplace=True)
df.head()
Out[20]:
In [21]:
mask = (df.index > 1250.0)
df = df[mask]
In [22]:
hdf = pd.read_csv("hetroscedastic.txt",names=['yerr', 'junk'], skiprows=1)
In [23]:
from etalon import *
np.random.seed(123456)
# Introduce the Real data
x = df.index.values
y_True = df.VG03_pos2.values/df.VG06.values
#Truncate to just 1200 nm
wl_mask = (x > 1200.0)
x = x[wl_mask]
y_True = y_True[wl_mask]
N = len(x)
# Define T_DSP for the model
T_DSP = T_gap_Si(x, 0.0)
#Introduce some noise with both measurement uncertainties
# and non-trivial correlated errors.
#yerr = 0.0002*np.ones(N)
yerr = hdf.yerr
iid_cov = np.diag(yerr ** 2)
a_True = 0.002
s_True = 100.0
y = y_True
In [24]:
plt.errorbar(x, y, yerr=yerr, fmt=".k", capsize=0)
plt.xlabel('$\lambda$');
plt.ylabel('$T_{gap}$');
plt.title('Measured Si Transmission');
plt.ylim(0.9, 1.01);
The off-diagonal terms are characterized by parameters $a$ and $s$.
The term $a$ controls the strength of the interaction. The term $s$ controls the range of correlation.
First we define a natural log likelihood function. You hand it a $d$, a $f$, an $a$, and an $s$ and it gives you the likelihood of the data: $\ln{p}=\ln{p(y|d,f,a,s)}$
Note: x is not passed to this function. The way Python works, x, y, and iid_cov are searched for in the local namespace, then global namespace.
In [25]:
n1 = sellmeier_Si(x)
In [26]:
def lnlike(d, f, lna, lns):
a, s = np.exp(lna), np.exp(lns)
off_diag_terms = a**2 * np.exp(-0.5 * (x[:, None] - x[None, :])**2 / s**2)
C = iid_cov + off_diag_terms
sgn, logdet = np.linalg.slogdet(C)
if sgn <= 0:
return -np.inf
r = y - T_gap_Si_withFF_fast(x, d, f, n1)/T_DSP
return -0.5 * (np.dot(r, np.linalg.solve(C, r)) + logdet)
In [27]:
# Apply a uniform prior over some range.
# Shape params are uniform in log space
def lnprior(d, f, lna, lns):
if not (3800 < d < 4300 and 0 < f < 0.5 and -12 < lna < -2 and 0 < lns < 8):
return -np.inf
return 0.0
In [28]:
def lnprob(p):
lp = lnprior(*p)
if not np.isfinite(lp):
return -np.inf
return lp + lnlike(*p)
In [29]:
ndim, nwalkers = 4, 32
d_Guess = 4120.0
f_Guess = 0.045
p0 = np.array([d_Guess, f_Guess, np.log(a_True), np.log(s_True)]) # Decent guess
pos = [p0 + 1.0e-2*p0 * np.random.randn(ndim) for i in range(nwalkers)]
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob)
In [30]:
@pickle_results('SiGaps_03_data-sampler.pkl')
def hammer_time(ndim, nwalkers, d_Guess, f_Guess, a_Guess, s_Guess, nburnins, ntrials):
# Initialize the walkers
p0 = np.array([d_Guess, f_Guess, np.log(a_Guess), np.log(s_Guess)])
pos = [p0 + 1.0e-2*p0 * np.random.randn(ndim) for i in range(nwalkers)]
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob)
pos, lp, state = sampler.run_mcmc(pos, nburnins)
sampler.reset()
pos, lp, state = sampler.run_mcmc(pos, ntrials)
return sampler
In [31]:
np.random.seed(78704)
ndim, nwalkers = 4, 32
d_Guess = 4120.0
f_Guess = 0.045
a_Guess = 0.0016
s_Guess = 100.0
nburnins = 200
ntrials = 600
In [32]:
sampler = hammer_time(ndim, nwalkers, d_Guess, f_Guess, a_Guess, s_Guess, nburnins, ntrials)
In [33]:
fig, axes = plt.subplots(4, 1, figsize=(5, 6), sharex=True)
fig.subplots_adjust(left=0.1, bottom=0.1, right=0.96, top=0.98,
wspace=0.0, hspace=0.05)
[a.plot(np.arange(chain.shape[1]), chain[:, :, i].T, "k", alpha=0.5)
for i, a in enumerate(axes)]
[a.set_ylabel("${0}$".format(l)) for a, l in zip(axes, ["d", "f", "\ln a", "\ln s"])]
axes[-1].set_xlim(0, chain.shape[1])
axes[-1].set_xlabel("iteration");
In [ ]:
samples = sampler.flatchain
samples_lin = copy(samples)
samples_lin[:, 2:] = np.exp(samples_lin[:, 2:])
In [ ]:
fig = triangle.corner(samples_lin,
labels=map("${0}$".format, ["d", "f", "a", "s"]),
quantiles=[0.16, 0.84])
In [ ]:
off_diag_terms = 0.0018**2 * np.exp(-0.5 * (x[:, None] - x[None, :])**2 / 12.0**2)
C = iid_cov + off_diag_terms
fit = T_gap_Si_withFF_fast(x, 4093.5, 0.04572, n1)/T_DSP
vec = np.random.multivariate_normal(fit, C)
In [ ]:
for d, f, a, s in samples_lin[np.random.randint(len(samples), size=3)]:
#plt.plot(x, T_gap_Si_withFF_fast(x, d, f, n1)/T_DSP, color="r", alpha=0.03)
off_diag_terms = a**2 * np.exp(-0.5 * (x[:, None] - x[None, :])**2 / s**2)
C = iid_cov + off_diag_terms
fit = T_gap_Si_withFF_fast(x, 4093.5, 0.04572, n1)/T_DSP
vec = np.random.multivariate_normal(fit, C)
plt.plot(x, vec, color="g", alpha=0.3)
plt.errorbar(x, y, yerr, fmt=".k")
plt.xlabel('$\lambda$');
plt.ylabel('$T_{gap}$');
plt.title('Observed Si Transmission');
plt.ylim(0.9, 1.01);
In [ ]:
d_mcmc, f_mcmc, a_mcmc, s_mcmc = map(lambda v: (v[1], v[2]-v[1], v[1]-v[0]),
zip(*np.percentile(samples_lin, [16, 50, 84],
axis=0)))
In [ ]:
d_mcmc, f_mcmc, a_mcmc, s_mcmc
In [ ]: