In [1]:
from __future__ import division, print_function
In [2]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
In [21]:
import george
from george import kernels
import pandas as pd
from astropy.stats import median_absolute_deviation as MAD
import glob
import emcee
import tqdm
In [22]:
# fn = '../data/215900519P-ep20160806.csv'
fns = glob.glob('../data/*.csv')
print(np.shape(fns))
fn = fns[4]
df = pd.read_csv(fn, skiprows=39)
df['time'] = df.t - df.t[0]
df['flux'] = df.fdt_t_roll_2D/ np.median(df.fdt_t_roll_2D) - 1.0
df['ferr'] = np.ones_like(df.t) * MAD(df.flux) /20
print(np.shape(df.time), np.shape(df.flux), np.shape(df.ferr))
In [23]:
plt.scatter(df.time,df.flux)
# plt.xlim(0,5)
plt.ylim(-0.03,0.03)
Out[23]:
In [24]:
y = df.loc[(df.time < 15) * (np.abs(df.flux)<0.02),'flux']
yerr = df.loc[(df.time < 15) * (np.abs(df.flux)<0.02),'ferr']
t = df.loc[(df.time < 15) * (np.abs(df.flux)<0.02),'time']
print(np.shape(t)), print(np.shape(yerr)), print(np.shape(t))
Out[24]:
In [25]:
A = 0.001
l = 10.
G = 2.0
sigma = 0.001
P = 2.0
kernel= A * kernels.ExpSquaredKernel(l) * kernels.ExpSine2Kernel(G, P) + kernels.WhiteKernel(sigma)
# kernel= A * kernels.ExpSine2Kernel(G, P) + kernels.WhiteKernel(sigma)
In [26]:
gp = george.GP(kernel, )#solver=george.HODLRSolver)
gp.compute(t, yerr=yerr)
print(gp.lnlikelihood(y))
print(gp.grad_lnlikelihood(y))
In [27]:
import scipy.optimize as op
# Define the objective function (negative log-likelihood in this case).
def nll(p, yval=y):
# Update the kernel parameters and compute the likelihood.
gp.kernel[:] = p
ll = gp.lnlikelihood(yval, quiet=True)
# The scipy optimizer doesn't play well with infinities.
return -ll if np.isfinite(ll) else 1e25
# And the gradient of the objective function.
def grad_nll(p, yval=y):
# Update the kernel parameters and compute the likelihood.
gp.kernel[:] = p
return -gp.grad_lnlikelihood(yval, quiet=True)
# You need to compute the GP once before starting the optimization.
gp.compute(t, yerr=yerr)
# Print the initial ln-likelihood.
print(gp.lnlikelihood(y))
# Run the optimization routine.
p0 = gp.kernel.vector
results = op.minimize(nll, p0, jac=grad_nll)
# Update the kernel and print the final log-likelihood.
gp.kernel[:] = results.x
print(gp.lnlikelihood(y))
In [28]:
print(results)
In [29]:
x = np.linspace(0, 15, 500)
mu, cov = gp.predict(y, x)
std = np.sqrt(np.diag(cov))
In [30]:
mup, covp = gp.predict(y, t + (t*0.0001)) # for some reason trying to predict at exactly the same time is a problem
stdp = np.sqrt(np.diag(covp))
flares = y > mup + (stdp*1.5)
tflares = np.where(flares,t,np.nan)
yflares = np.where(flares,y,np.nan)
In [31]:
fig, ax = plt.subplots(1,1,figsize=[11,6], )
ax.errorbar(t, y, yerr=yerr, fmt=".k", capsize=0,alpha=0.7)
ax.plot(tflares, yflares, 'o-')
ax.fill_between(x,mu-std,mu+std, color="#ff7f0e", alpha=0.3,
edgecolor="none")
ax.plot(x,mu,color="#ff7f0e", alpha=0.7)
ax.set_xlim(8,15)
ax.set_ylabel('Relative brightness', fontsize=17)
ax.minorticks_on()
In [32]:
gp = george.GP(kernel, )#solver=george.HODLRSolver)
gp.compute(t[~flares], yerr=yerr[~flares])
print(gp.lnlikelihood(y[~flares]))
print(gp.grad_lnlikelihood(y[~flares]))
# You need to compute the GP once before starting the optimization.
gp.compute(t[~flares])
# Print the initial ln-likelihood.
print(gp.lnlikelihood(y[~flares]))
# Run the optimization routine.
p0 = gp.kernel.vector
results = op.minimize(nll, p0, jac=grad_nll,args=y[~flares])
# Update the kernel and print the final log-likelihood.
gp.kernel[:] = results.x
print(gp.lnlikelihood(y[~flares]))
In [33]:
x = np.linspace(0, 15, 500)
mu, cov = gp.predict(y[~flares], x)
std = np.sqrt(np.diag(cov))
In [34]:
mup, covp = gp.predict(y[~flares], t + (t*0.0001)) # for some reason trying to predict at exactly the same time is a problem
stdp = np.sqrt(np.diag(covp))
flares = y > mup + (stdp*2.5)
tflares = np.where(flares,t,np.nan)
yflares = np.where(flares,y,np.nan)
In [35]:
fig, ax = plt.subplots(1,1,figsize=[11,6], )
ax.errorbar(t, y, yerr=0, fmt=".k", capsize=0,alpha=0.7)
ax.plot(tflares, yflares, 'o')
ax.fill_between(x,mu-std*2,mu+std*2, color="#ff7f0e", alpha=0.3,
edgecolor="none")
ax.plot(x,mu,color="#ff7f0e", alpha=0.7)
ax.set_xlim(8,15)
ax.set_ylabel('Relative brightness', fontsize=17)
ax.minorticks_on()
In [36]:
(np.e**gp.kernel[:])
Out[36]:
In [37]:
gp.kernel.vector
Out[37]:
In [38]:
nwalkers, ndim = 16, len(gp.kernel.vector)
In [85]:
def angus_kernel(theta):
"""
use the kernel that Ruth Angus uses. Be sure to cite her
"""
theta = np.exp(theta)
A = theta[0]
l = theta[1]
G = theta[2]
sigma = theta[4]
P = theta[3]
kernel = (A * kernels.ExpSquaredKernel(l) *
kernels.ExpSine2Kernel(G, P) +
kernels.WhiteKernel(sigma)
)
return kernel
def lnprob(p, time, y, yerr):
# Trivial improper prior: uniform in the log.
# from DFM george user guide
if np.any((-p < -20) + (p > 10)):
return -np.inf
lnprior = 0.0
kernel = angus_kernel(p)
gp = george.GP(kernel)
gp.compute(time, yerr)
return lnprior + gp.lnlikelihood(y, quiet=True)
In [86]:
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob,
args=[
t,
y,
yerr,
])
In [87]:
p0 = (np.repeat(gp.kernel.vector, nwalkers) *
np.random.random(size=ndim * nwalkers))
p0 = p0.reshape(ndim, nwalkers).T
In [88]:
p0, _, _ = sampler.run_mcmc(p0, 1)
In [128]:
samples = sampler.run_mcmc(p0, 2000)
In [151]:
plt.hist(np.e**sampler.chain[:,1500:,3].flatten())
Out[151]:
In [158]:
Out[158]:
In [ ]:
In [ ]:
In [ ]:
In [ ]: