This IPython Notebook is for performing a fit and generating a figure of the spectrum of sample VG09-12, in the mesh region with 49+/-6 nm gap. This version is modified to fit for two gaps, updated with new data from February 2015.

The filename of the figure is [TBD].pdf.

Author: Michael Gully-Santiago, gully@astro.as.utexas.edu

Date: February 16, 2015


In [1]:
%pylab inline
import emcee
import triangle
import pandas as pd
import seaborn as sns
from astroML.decorators import pickle_results


Populating the interactive namespace from numpy and matplotlib

In [2]:
sns.set_context("paper", font_scale=2.0, rc={"lines.linewidth": 2.5})
sns.set(style="ticks")

In [3]:
from etalon import *
np.random.seed(78704)

In [4]:
df = pd.read_csv('../data/VG09_12_gap_20150206.csv', index_col=0)
df.head()


Out[4]:
VG09-12_1_0 VG09-12_1_2 VG09-12_1_4 VG09-12_1_6 VG09-12_1_8 VG09-12_1_10 VG09-12_1_12 VG09-12_1_14 VG09-12_1_16 VG09-12_1_18 ... VG09-12_1_32 VG09-12_1_34 VG09-12_1_36 VG09-12_1_38 VG09-12_1_40 VG09-12_1_42 VG09-12_1_44 VG09-12_1_46 VG09-12_1_48 VG09-12_1_50
wavelength
1775 1.880152 1.641429 0.976230 0.964145 0.976460 0.979872 0.987270 0.986529 1.070819 1.002500 ... 1.002960 1.002864 1.003209 1.003811 0.965806 0.918209 1.876167 1.881628 1.883350 1.883947
1765 1.879758 1.642331 0.974838 0.964183 0.976321 0.978982 0.986934 0.985437 1.070772 1.002522 ... 1.002708 1.002673 1.003248 1.003872 0.966022 0.917161 1.875934 1.881610 1.883944 1.884051
1755 1.880479 1.643144 0.973691 0.964016 0.975474 0.978551 0.986437 0.984843 1.070936 1.002670 ... 1.002381 1.002297 1.003009 1.003320 0.965691 0.917684 1.876912 1.882451 1.884022 1.884329
1745 1.880220 1.643559 0.973281 0.963117 0.974474 0.977558 0.986889 0.984795 1.070851 1.002239 ... 1.002931 1.002059 1.002659 1.003579 0.965511 0.916989 1.876531 1.882771 1.884626 1.884169
1735 1.880693 1.643988 0.973618 0.963056 0.973207 0.976588 0.987083 0.984501 1.070377 1.001896 ... 1.002166 1.002575 1.003052 1.003299 0.965231 0.917656 1.876819 1.882872 1.884520 1.884540

5 rows × 26 columns


In [5]:
# Introduce the Real data, decimate the data.
x = df.index.values
N = len(x)
# Define T_DSP for the model
T_DSP = T_gap_Si(x, 0.0)
n1 = sellmeier_Si(x)

# Define uncertainty
yerr = 0.0002*np.ones(N)
yerr[(x > 1350) & (x < 1420)] = 0.0005 #higher noise in this region
iid_cov = np.diag(yerr ** 2)

# Select the spectrum of interest
# Normalize the spectrum by measured DSP Si wafer.
y = df['VG09-12_1_26'].values

Define the likelihood.


In [6]:
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)

Define the prior.


In [7]:
def lnprior(d, f, lna, lns):
    if not (0 < d < 100 and 0.0 < f < 1.0 and -12 < lna < -2 and 0 < lns < 10):
        return -np.inf
    return 0.0

Combine likelihood and prior to obtain the posterior.


In [8]:
def lnprob(p):
    lp = lnprior(*p)
    if not np.isfinite(lp):
        return -np.inf
    return lp + lnlike(*p)

Set up emcee.


In [9]:
@pickle_results('SiGaps_18_0gap.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

Set up the initial conditions.


In [10]:
np.random.seed(78704)
ndim, nwalkers = 4, 32
d_Guess = 15.0
f_Guess = 0.95
a_Guess = 0.0016
s_Guess = 500.0
nburnins = 300
ntrials = 9000

Run the burn-in phase. Run the full MCMC. Pickle the results.


In [11]:
sampler = hammer_time(ndim, nwalkers, d_Guess, f_Guess, a_Guess, s_Guess, nburnins, ntrials)


@pickle_results: using precomputed results from 'SiGaps_18_0gap.pkl'

Linearize $a$ and $s$ for easy inspection of the values.


In [12]:
chain = sampler.chain
samples_lin = copy(sampler.flatchain)
samples_lin[:, 2:] = np.exp(samples_lin[:, 2:])

In [13]:
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");


Make a triangle corner plot.


In [14]:
fig = triangle.corner(samples_lin, 
                      labels=map("${0}$".format, ["d", "f", "a", "s"]), 
                      quantiles=[0.16, 0.84])


Quantiles:
[(0.16, 3.9775785423732453), (0.84, 31.600191065182045)]
Quantiles:
[(0.16, 0.03449491294438941), (0.84, 0.7088700630482242)]
Quantiles:
[(0.16, 0.0010558971716583014), (0.84, 0.0096825717902464681)]
Quantiles:
[(0.16, 311.06441219384976), (0.84, 4062.0138060398622)]

In [15]:
fig = triangle.corner(samples_lin[:,0:2], 
                      labels=map("${0}$".format, ["d", "f"]), 
                      quantiles=[0.16, 0.84])
plt.savefig("VG0912_0gap_corner.pdf")


Quantiles:
[(0.16, 3.9775785423732453), (0.84, 31.600191065182045)]
Quantiles:
[(0.16, 0.03449491294438941), (0.84, 0.7088700630482242)]
/Users/gully/anaconda/lib/python2.7/site-packages/matplotlib/backends/backend_pdf.py:2264: FutureWarning: comparison to `None` will result in an elementwise object comparison in the future.
  different = bool(ours != theirs)

In [16]:
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)))
d_mcmc, f_mcmc, a_mcmc, s_mcmc


Out[16]:
((10.677500418273841, 20.922690646908205, 6.6999218759005954),
 (0.24973230076548836, 0.45913776228273584, 0.21523738782109894),
 (0.0026829125467427973, 0.0069996592435036704, 0.0016270153750844959),
 (956.47471519625674, 3105.5390908436057, 645.41030300240698))

In [17]:
print "{:.0f}^{{+{:.0f}}}_{{-{:.0f}}}".format(*d_mcmc)
print "{:.3f}^{{+{:.3f}}}_{{-{:.3f}}}".format(*f_mcmc)


11^{+21}_{-7}
0.250^{+0.459}_{-0.215}

In [26]:
plt.figure(figsize=(4,3))

for d, f, a, s in samples_lin[np.random.randint(len(samples_lin), size=60)]:
    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, d, f, n1)/T_DSP
    vec = np.random.multivariate_normal(fit, C)
    plt.plot(x, vec,"-b", alpha=0.06)

plt.step(x, y,color="k", label='Measurement')
plt.errorbar(x, y,yerr=yerr, color="k")
fit = T_gap_Si_withFF_fast(x, 14, 1.0, n1)/T_DSP
fit_label = 'Model with $d={:.0f}$ nm, $f={:.2f}$'.format(14, 1.0)
plt.plot(x, fit, '--', color=sns.xkcd_rgb["pale red"], alpha=1.0, label=fit_label)
plt.plot([-10, -9], [-10, -9],"-b", alpha=0.45, label='Draws from GP')
plt.plot([0, 5000], [1.0, 1.0], '-.k', alpha=0.5)
plt.fill_between([1200, 1250], 2.0, 0.0, hatch='\\', alpha=0.4, color='k', label='Si absorption cutoff')

plt.xlabel('$\lambda$ (nm)');
plt.ylabel('$T_{gap}$');
plt.xlim(1200, 1900);
plt.ylim(0.96, 1.019);
plt.legend(loc='lower right')
plt.savefig("VG0912_0gap.pdf",  bbox_inches='tight')



In [19]:
y.mean()


Out[19]:
0.99980907280714904

In [20]:
y.std()


Out[20]:
0.00084827408755901135

In [21]:
ids = samples_lin[:,1] > 0.8
ids.sum()/(len(ids)*1.0)
np.percentile(samples_lin[ids, 0], [99.5])


Out[21]:
array([ 13.23949635])

The VG09-12 off-mesh spectrum rules out gaps greater than 13 nm over $>80$\% of the measurement area.