In [1]:
import os
import numpy as np
import astropy
import astropy.units as u
print astropy.__version__
%matplotlib inline
import matplotlib.pyplot as plt
In [2]:
from JLA import PhotometryData
In [3]:
fname = '/Users/lisaleemcb/data/jla_light_curves/lc-03D1au.list'
In [4]:
jla = PhotometryData(fname)
In [5]:
lc = jla.photometryTable
lc_ext = jla.photometryTable
In [6]:
import sncosmo
In [7]:
bandPassList = ['u', 'g', 'r', 'i', 'z']
banddir = os.path.join(os.getenv('THROUGHPUTS_DIR'), 'megacam')
In [8]:
for band in bandPassList:
# setup sncosmo bandpasses
bandfname = banddir + '/' + band + 'Mega' + '.fil.txt'
print bandfname
numpyband = np.loadtxt(bandfname)
sncosmoband = sncosmo.Bandpass(wave=numpyband[:, 0],
trans=numpyband[:, 1],
wave_unit=u.nm,
name='SNLS_' + band)
sncosmo.registry.register(sncosmoband, force=True)
In [9]:
source = sncosmo.get_source('salt2-extended')
dust = sncosmo.CCM89Dust()
model = sncosmo.Model(source=source,
effects=[dust, dust],
effect_names=['host', 'mw'],
effect_frames=['rest', 'obs'])
params = {'z': 0.50308369405, 't0': 52909.7452198, 'x0': 1.13167971126e-05, 'x1': 1.2731910211, 'c': -0.0123529588226}
model.set(**params)
model_ext = sncosmo.Model(source=source,
effects=[dust, dust],
effect_names=['host', 'mw'],
effect_frames=['rest', 'obs'])
model_ext.set(**params)
model_ext.set(hostebv=0.05)
In [10]:
fig_jla = sncosmo.plot_lc(lc)
In [11]:
lc['modelflux'] = model.bandflux(time=lc['time'], band=lc['band'], zp=lc['zp'], zpsys=lc['zpsys'])
lc['deviation'] = np.random.normal(0,1, size=len(lc)) * lc['fluxerr']
lc_ext['modelflux'] = model_ext.bandflux(time=lc_ext['time'], band=lc_ext['band'], zp=lc_ext['zp'],
zpsys=lc_ext['zpsys'])
lc_ext['deviation'] = np.random.normal(0,1, size=len(lc)) * lc_ext['fluxerr']
In [12]:
lc['scatteredFlux'] = lc['modelflux'] + lc['deviation']
lc_ext['scatteredFlux'] = lc_ext['modelflux'] + lc_ext['deviation']
In [13]:
lc_ext[:5]
Out[13]:
In [14]:
# rename the columns so sncosmo pick the right ones for analysis
lc.rename_column('flux' , 'jla_flux')
lc.rename_column('scatteredFlux', 'flux')
lc_ext.rename_column('flux', 'jla_flux')
lc_ext.rename_column('scatteredFlux', 'flux')
lc[:5]
Out[14]:
In [73]:
fig_modelflux = sncosmo.plot_lc(lc, model=[model], pulls=False, show_model_params=False)
fig_modelflux.savefig('lc.png')
In [16]:
import triangle
In [17]:
mcmc_res_noext, mcmc_model_noext = sncosmo.mcmc_lc(lc, model,
vparam_names=['t0', 'x0', 'x1', 'c'],
bounds={'c':(-0.3, 0.3), 'x1':(-3.0, 3.0)}, minsnr=3.0)
In [18]:
mcmc_res, mcmc_model = sncosmo.mcmc_lc(lc, model, vparam_names=['t0', 'x0', 'x1', 'c', 'hostebv'],
bounds={'c':(-0.3, 0.3), 'x1':(-3.0, 3.0)}, minsnr=3.0)
In [19]:
mcmc_res_ext, mcmc_model_ext = sncosmo.mcmc_lc(lc_ext, model_ext,
vparam_names=['t0', 'x0', 'x1', 'c', 'hostebv'],
bounds={'c':(-0.3, 0.3), 'x1':(-3.0, 3.0)}, minsnr=3.0)
In [20]:
mcmc = mcmc_res.vparam_names
mcmc_noext = mcmc_res_noext.vparam_names
mcmc_ext = mcmc_res_ext.vparam_names
mcmc_noext_samples = mcmc_res_noext.samples
mcmc_samples = mcmc_res.samples
mcmc_ext_samples = mcmc_res_ext.samples
mcmc
Out[20]:
In [35]:
alpha = 0.11
beta = -3.1
mB = source.peakmag(band='bessellB', magsys='ab')
print mB
noext_dict = {}
mcmc_dict = {}
ext_dict = {}
In [37]:
noext_dict['x0'] = mcmc_noext_samples[:,1]
noext_dict['x1'] = mcmc_noext_samples[:,2]
noext_dict['c'] = mcmc_noext_samples[:,3]
noext_dict['mu'] = mB* + alpha * noext_dict['x1'] + beta * noext_dict['c']
In [57]:
mcmc_dict['x0'] = mcmc_samples[:,1]
mcmc_dict['x1'] = mcmc_samples[:,2]
mcmc_dict['c'] = mcmc_samples[:,3]
mcmc_dict['hostebv'] = mcmc_samples[:,4]
mcmc_dict['mu'] = mB* + alpha * mcmc_dict['x1'] + beta * mcmc_dict['c']
In [58]:
ext_dict['x0'] = mcmc_ext_samples[:,1]
ext_dict['x1'] = mcmc_ext_samples[:,2]
ext_dict['c'] = mcmc_ext_samples[:,3]
ext_dict['hostebv'] = mcmc_ext_samples[:,4]
ext_dict['mu'] = mB* + alpha * ext_dict['x1'] + beta * ext_dict['c']
In [61]:
noext_arr = []
mcmc_arr =[]
ext_arr = []
In [62]:
noext_arr.append(noext_dict['c'])
noext_arr.append(noext_dict['mu'])
mcmc_arr.append(mcmc_dict['c'])
mcmc_arr.append(mcmc_dict['hostebv'])
mcmc_arr.append(mcmc_dict['mu'])
ext_arr.append(ext_dict['c'])
ext_arr.append(ext_dict['hostebv'])
ext_arr.append(ext_dict['mu'])
In [49]:
np.shape(np.transpose(noext_arr))
Out[49]:
In [68]:
figure = triangle.corner(np.transpose(noext_arr), labels=['c', r'$\mu$'],
truths=[model.get('c'), mB* + alpha * model.get('x1') + beta * model.get('c')],
range=2*[0.9999], show_titles=True, title_args={"fontsize": 12},
hist_kwargs={'normed':True})
figure.gca().annotate(" ", xy=(0.5, 1.0), xycoords="figure fraction",
xytext=(0, -5), textcoords="offset points",
ha="center", va="top")
figure.savefig('noext.png')
In [69]:
figure = triangle.corner(np.transpose(mcmc_arr), labels=['c','hostebv', r'$\mu$'],
truths=[model.get('c'), model.get('hostebv'),
mB* + alpha * model.get('x1') + beta * model.get('c')],
range=3*[0.9999], show_titles=True, title_args={"fontsize": 12},
hist_kwargs={'normed':True})
figure.gca().annotate("mcmc sampling host_ebv=0", xy=(0.5, 1.0), xycoords="figure fraction",
xytext=(0, -5), textcoords="offset points",
ha="center", va="top")
figure.savefig('mcmc.png')
In [70]:
figure = triangle.corner(np.transpose(ext_arr), labels=['c','hostebv', r'$\mu$'],
truths=[model_ext.get('c'), model_ext.get('hostebv'),
mB* + alpha * model_ext.get('x1') + beta * model_ext.get('c')],
range=3*[0.9999], show_titles=True, title_args={"fontsize": 12},
hist_kwargs={'normed':True})
figure.gca().annotate("mcmc sampling host_ebv=.05", xy=(0.5, 1.0), xycoords="figure fraction",
xytext=(0, -5), textcoords="offset points",
ha="center", va="top")
figure.savefig('ext.png')
In [42]:
mcmc_samples[:,3].min(), mcmc_samples[:,3].max()
Out[42]:
In [129]:
import plotly.plotly as py
import plotly.graph_objs as go
import numpy as np
t0 = mcmc_samples[:,0]
x0 = mcmc_samples[:,1]
x1 = mcmc_samples[:,2]
c = mcmc_samples[:,3]
hostebv = mcmc_samples[:,4]
trace1 = go.Scatter(
x=x, y=y, mode='markers', name='points',
marker=dict(color='gray', size=2, opacity=0.4)
)
trace2 = go.Histogram2dcontour(
x=x, y=y, name='density', ncontours=10,
colorscale=[[0, 'rgb(255,255,255)'],
[.25, 'darkorange'],
[0.5, 'yellow'],
[0.75, 'green'],
[1,'blue']],
reversescale=False, showscale=True
)
histt0 = go.Histogram(
x=t0, name='x density', opacity=0.5,
marker=dict(color='blue'),
xaxis='x', yaxis='y'
)
histx0 = go.Histogram(
x=x0, name='y density', opacity=0.5,
marker=dict(color='blue'),
xaxis='x2', yaxis='y3'
)
histx1 = go.Histogram(
x=x1, name='x density', opacity=0.5,
marker=dict(color='blue'),
xaxis='x3'
)
histc = go.Histogram(
x=c, name='y density', opacity=0.5,
marker=dict(color='blue'),
xaxis='x4'
)
data = [histt0, histx0]
layout = go.Layout(
showlegend=False,
autosize=False,
width=600,
height=550,
xaxis=dict(
domain=[0, 0.2],
showgrid=False,
zeroline=False
),
yaxis=dict(
domain=[.8, 1],
showgrid=False,
zeroline=False
),
margin=dict(
t=50
),
hovermode='closest',
bargap=0,
xaxis2=dict(
domain=[0.2, 1],
showgrid=False,
zeroline=False
),
yaxis2=dict(
domain=[0, 0.2],
showgrid=False,
zeroline=False
),
xaxis3=dict(
domain=[0.2, 1],
showgrid=False,
zeroline=False
),
yaxis3=dict(
domain=[0.2, 1],
showgrid=False,
zeroline=False
)
)
iplot(go.Figure(data=data, layout=layout))
In [ ]: