WASP-80b broadband analysis

3a. Gaussian process hyperparameter estimation I

Hannu Parviainen, Instituto de Astrofísica de Canarias


This notebook works as an appendix to Parviainen et al., Ground based transmission spectroscopy of WASP-80b (2017). The paper covers two analyses: a broadband analysis using three previously published datasets, and a transmission spectroscopy analysis using two GTC-observed spectroscopic time series, and this notebook covers a part of the broadband analysis.

Last (significant) revision: 11.08.2017


Here we estimate the GP hyperparameters (HPs) for the light curves in T13 and M14 datasets. The GP for these two datasets uses a simple kernel with time as the only input parameter.


In [1]:
%pylab inline
%run __init__.py


Populating the interactive namespace from numpy and matplotlib

In [2]:
from exotk.utils.misc import fold
from src.extcore import *

Compute residuals

The GP hyperparameters are fit to the residuals from the ckwn white-noise analysis.


In [3]:
lpf = LPFTM()
pv0 = pd.read_hdf(RFILE_EXT, 'ckwn/fc').median().values

fluxes_m = lpf.compute_transit(pv0)
residuals = [fo-fm for fo,fm in zip(lpf.fluxes, fluxes_m)]
gps = [GPTime(time, res) for time,res in zip(lpf.times, residuals)]
hps = []

Plot the light curves


In [4]:
phases = list(map(lambda t: fold(t, P, TC, 0.5)-0.5, lpf.times))
fig,axs = subplots(4,3, figsize=(14,14),sharey=True, sharex=True)
for iax,ilc in enumerate(lpf.lcorder):
    a = axs.flat[iax]
    a.plot(phases[ilc], lpf.fluxes[ilc],'.', alpha=0.5)
    a.plot(phases[ilc], fluxes_m[ilc],'k')
    a.plot(phases[ilc], lpf.fluxes[ilc]-fluxes_m[ilc]+0.95,'.', alpha=0.5)
    a.text(0.5, 0.95, lpf.passbands[ilc], ha='center', va='top', size=12, transform=a.transAxes)
setp(axs, ylim=(0.94,1.01), xlim=(-0.035,0.035))
fig.tight_layout()
axs.flat[-1].set_visible(False)


Fit the Hyperparameters and plot the GP mean with the data


In [5]:
hps = []
for gp in tqdm(gps, desc='Optimising GP hyperparameters'):
    gp.fit()
    hps.append(gp.hp)


Optimising GP hyperparameters: 100%|██████████| 11/11 [00:05<00:00,  3.01it/s]

In [6]:
fig,axs = subplots(4,3, figsize=(14,10),sharey=True, sharex=True)
for iax,ilc in enumerate(lpf.lcorder):
    axs.flat[iax].plot(phases[ilc], gps[ilc].flux, '.', alpha=0.5)
    gps[ilc].compute(hps[ilc])
    pr = gps[ilc].predict()
    axs.flat[iax].plot(phases[ilc], pr, 'k')
setp(axs, ylim=(-0.015,.015), xlim=(-0.04,0.04))
fig.tight_layout()
axs.flat[-1].set_visible(False)



In [7]:
fig,axs = subplots(4,3, figsize=(14,10),sharey=True, sharex=True)
for iax,ilc in enumerate(lpf.lcorder):
    axs.flat[iax].plot(phases[ilc], lpf.fluxes[ilc], '.', alpha=0.5)
    gps[ilc].compute(hps[ilc])
    pr = gps[ilc].predict()
    axs.flat[iax].plot(phases[ilc], fluxes_m[ilc]+pr, 'k')
setp(axs, ylim=(0.955,1.015), xlim=(-0.04,0.04))
fig.tight_layout()
axs.flat[-1].set_visible(False)


Create a Pandas dataframe and save the hyperparameters


In [8]:
with pd.HDFStore(DFILE_EXT) as f:
    ntr = [k[3:] for k in f.keys() if 'lc/triaud' in k]
    nma = [k[3:] for k in f.keys() if 'lc/mancini' in k]

df = pd.DataFrame(hps, columns=gp.names, index=lpf.passbands)
df['lc_name'] = ntr+nma
df


Out[8]:
ln_wn_var ln_output_var ln_input_scale lc_name
r -14.061009 -15.602592 -9.482634 /triaud2013/r/eulercam_1
z -11.538013 -14.091494 -11.006471 /triaud2013/z/trappist_1
z -11.326154 -15.159640 -6.415521 /triaud2013/z/trappist_2
H -11.992534 -12.533673 -10.651665 /mancini2014/H/grond
J -11.255762 -13.003909 -8.117159 /mancini2014/J/grond
K -10.953819 -11.730327 -9.770935 /mancini2014/K/grond
g -14.402214 -16.615069 -5.768119 /mancini2014/g/grond
i -14.380497 -30.701713 -1.322828 /mancini2014/i/dfosc
i -14.380497 -30.701713 -1.322828 /mancini2014/i/grond
r -15.508300 -17.281088 -9.295438 /mancini2014/r/grond
z -14.729695 -14.270637 -11.744267 /mancini2014/z/grond

In [9]:
df.ix[:3].to_hdf(RFILE_EXT, 'gphp/triaud2013')
df.ix[3:].to_hdf(RFILE_EXT, 'gphp/mancini2014')