In [1]:
%load_ext autoreload
In [2]:
import numpy as np
import pandas as pd
%matplotlib inline
import matplotlib.pyplot as plt
import os
In [3]:
#from OpSimSummary import summarize_opsim as oss
import sncosmo
In [4]:
#import cov_utils as cutils
In [5]:
# sncosmo Bandpasses required for fitting
throughputsdir = os.getenv('THROUGHPUTS_DIR')
from astropy.units import Unit
bandPassList = ['u', 'g', 'r', 'i', 'z', 'y']
banddir = os.path.join(os.getenv('THROUGHPUTS_DIR'), 'baseline')
for band in bandPassList:
# setup sncosmo bandpasses
bandfname = banddir + "/total_" + band + '.dat'
# register the LSST bands to the SNCosmo registry
# Not needed for LSST, but useful to compare independent codes
# Usually the next two lines can be merged,
# but there is an astropy bug currently which affects only OSX.
numpyband = np.loadtxt(bandfname)
print band
sncosmoband = sncosmo.Bandpass(wave=numpyband[:, 0],
trans=numpyband[:, 1],
wave_unit=Unit('nm'),
name=band)
sncosmo.registry.register(sncosmoband, force=True)
Read in the light curve
In [53]:
lc = sncosmo.read_lc('example_data/lsst_median.dat')
In [55]:
pd.DataFrame(np.asarray(lc)).head()
Out[55]:
Setup the SNCosmo model
In [58]:
model = sncosmo.Model(source='salt2-extended')
model.set(**{'z':0.5, 't0':49540, 'x0':1.0e-5})
print model
Plot the light curve
In [59]:
sncosmo.plot_lc(lc, model=model, color='k', pulls=False)
Out[59]:
Obtain the mcmc_output
In [14]:
mcmc_out = sncosmo.mcmc_lc(lc, model=model,
vparam_names=['t0', 'x0', 'x1', 'c'],
bounds={'c':(-1., 1.), 'x1':(-5.0, 5.0)},
minsnr=0., modelcov=False)
In [60]:
import analyzelcFits as anf
In [61]:
t = anf.ResChar.fromSNCosmoRes(mcmc_out)
Now we can look at the MCMC samples as they were
In [64]:
t.samples.head()
Out[64]:
We can convert to 'SAlt samples' which gives us a few more parameters In calculating mu, values of alpha, beta, MDelta are assumed
In [65]:
samples = t.salt_samples()
In [66]:
samples.head()
Out[66]:
Though we put in MDelta, this is currently 0, so we will remove it
In [67]:
x = samples[['t0', 'x0', 'x1', 'c', 'mB','mu']]
In [68]:
x.head()
Out[68]:
In [69]:
x.values
Out[69]:
In [70]:
from corner import corner
In [71]:
corner(x.values, labels=x.columns)
Out[71]:
We can see the error on mu for these values of alpha and beta
In [72]:
t.salt_samples().mu.std()
Out[72]:
This is the covariance matrix that was read in as output of mcmc
In [73]:
t._covariance
Out[73]:
We display it in a more easy to read format
In [74]:
t.covariance
Out[74]:
An advantage is that we can now get a covariance element very easily either by index or by parameter
In [78]:
t.covariance.ix['x0', 'x0']
Out[78]:
In [79]:
t.covariance.ix['x0', 'x1']
Out[79]:
In [80]:
t.covariance.iloc[1,2]
Out[80]:
For submatrices
In [81]:
t.covariance.ix[['x0', 'x1'], ['x0', 'x1']]
Out[81]:
In [82]:
import cov_utils as cu
In [83]:
cu.subcovariance(t.covariance, ['x0', 'c'])
Out[83]:
We can calculate the covariance of the parameters that go into the distance modulus formula using a linear approximation
In [75]:
m = t.salt_covariance_linear(); m
Out[75]:
We could have also used the mcmc samples to estimate the covariance without resorting to assumptions on linearity
In [76]:
ss = t.salt_samples()[['mB', 'x1', 'c']]
In [77]:
ss.cov()
Out[77]:
Which match reasonably well
In [84]:
import seaborn as sns
In [85]:
sns.set()
In [86]:
sns.jointplot(x='mu', y='c', data=samples, kind='kde')
Out[86]:
In [87]:
from scipy import stats
In [88]:
sns.distplot(samples.mu, hist=True, rug=False, hist_kws={'normed':1, 'histtype':'step',
'lw':1, 'alpha':1, 'color':'k'},
fit=stats.norm)
Out[88]:
In [89]:
del samples['MDelta']
In [90]:
samples.head()
Out[90]:
In [92]:
del samples['x0']
In [93]:
h = sns.PairGrid(samples);
h.map_upper(plt.scatter)
h.map_diag(sns.distplot, rug=False, hist_kws={'histtype':'step', 'alpha':1, 'lw':1})
h.map_lower(sns.kdeplot, legend=False, lw=3, n_levels=3)
In [ ]:
h.map(plt.scatter)
In [ ]: