In [107]:
import numpy as np
import pandas as pd
import os
import time
In [2]:
%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()
In [3]:
import sncosmo
import analyzeSN as ans
Download the simulations located at http://lsst.astro.washington.edu/simdata/SN_data/MINION/MINION_1016_SNANA_DDF_10yr_v2.tgz and un-tar them
tar -xzvf simulationDirFor me the resultant directory is at
/Users/rbiswas/data/LSST/SNANA_data/MINION_1016_10YR_DDF_v2
In [4]:
!tree /Users/rbiswas/data/LSST/SNANA_data/MINION_1016_10YR_DDF_v2/
In [5]:
simdata = ans.SNANASims.fromSNANAfileroot('LSST_Ia',
location='/Users/rbiswas/data/LSST/SNANA_data/MINION_1016_10YR_DDF_v2/',
coerce_inds2int=False)
In [6]:
simdata.headData.head()
Out[6]:
In [7]:
hd = simdata.headData.query('REDSHIFT_FINAL > 0.2 and REDSHIFT_FINAL < 0.25').head(5).copy()
In [8]:
hd
Out[8]:
In [9]:
lcobj = simdata.get_SNANA_photometry(snid='2341')
In [10]:
simdata.headData.ix['2341']
Out[10]:
In [11]:
from collections import OrderedDict as odict
In [12]:
hd['SIM_SALT2x0'] = hd.SIM_SALT2x0 * 10.0**(-0.4 * 0.27)
In [13]:
hd['params'] = hd[['SIM_PEAKMJD', 'SIM_SALT2x0', 'SIM_SALT2x1', 'SIM_SALT2c', 'REDSHIFT_FINAL', ]].rename(columns=dict(SIM_PEAKMJD='t0',
SIM_SALT2x0='x0',
SIM_SALT2x1='x1',
SIM_SALT2c='c',
SIM_MWEBV='ebv',
REDSHIFT_FINAL='z')).apply(odict, axis=1)
In [14]:
dust = sncosmo.CCM89Dust()
model = sncosmo.Model(source='salt2-extended', effects=[dust, dust],
effect_names=['host', 'mw'], effect_frames=['rest', 'obs'])
model.set(**hd.ix['2341'].params)
Here we compare the data to the model
In [15]:
_ = sncosmo.plot_lc(lcobj.snCosmoLC(), model=model)
First we will use a max likelihood method for estimating the model parameters
In [16]:
fitres = sncosmo.fit_lc(lcobj.snCosmoLC(), model=model, vparam_names=['t0', 'x0', 'x1', 'c'], modelcov=True)
In [17]:
fitresChar = ans.ResChar.fromSNCosmoRes(fitres)
In [18]:
fitresChar.salt_covariance_linear()
Out[18]:
In [19]:
fitresChar.mu_variance_linear()** 0.5
Out[19]:
In [20]:
fitres_z = sncosmo.fit_lc(lcobj.snCosmoLC(), model=model, vparam_names=['t0', 'x0', 'x1', 'c', 'z'], guess_z=True,
bounds=dict(z=(0.1, 0.35)),
modelcov=True)
In [21]:
fitresChar_z = ans.ResChar.fromSNCosmoRes(fitres_z)
In [22]:
fitresChar_z.mu_variance_linear()**0.5
Out[22]:
The maximum likelihood method above is not good if the likelihood are too far from the Gaussian. A method to sample posteriors is better. Such a method can be based on a MCMC run. We do that here by the following code
In [24]:
samples = sncosmo.mcmc_lc(lcobj.snCosmoLC(), model=model, vparam_names=['t0', 'x0', 'x1', 'c'], modelcov=True)
These are samples of these parameters :
In [50]:
ans.ResChar.fromSNCosmoRes(samples).samples.head()
Out[50]:
In [108]:
tstart = time.time()
samples_z_nomodelcov = sncosmo.mcmc_lc(lcobj.snCosmoLC(), model=model, vparam_names=['z', 't0', 'x0', 'x1', 'c'],
bounds=dict(z=(0.1, 0.5)), modelcov=False)
tend = time.time()
print(tend-tstart)
In [110]:
tstart = time.time()
samples_z = sncosmo.mcmc_lc(lcobj.snCosmoLC(), model=model, vparam_names=['z', 't0', 'x0', 'x1', 'c'],
bounds=dict(z=(0.2, 0.6)), modelcov=True)
tend = time.time()
print(tend-tstart)
In [52]:
ans.ResChar.fromSNCosmoRes(samples_z).samples.head()
Out[52]:
In [57]:
fig, ax = plt.subplots()
ans.ResChar.fromSNCosmoRes(samples_z).samples.z.hist(histtype='step', bins=20, lw=2., alpha=1., normed=1, ax=ax)
ax.axvline(0.22)
Out[57]:
In [101]:
import corner
The origin of the bimodal distribution on mu can be seen by exploring the scatter plots of the samples.
For the case where z was known, we can look at the histograms and scatter plots of pairs of variables:
In [100]:
sns.pairplot(ans.ResChar.fromSNCosmoRes(samples).salt_samples(), kind='scatter', diag_kind='kde')
Out[100]:
Consequently, the posterior distributions can be plotted and they are fine
In [102]:
corner.corner(ans.ResChar.fromSNCosmoRes(samples).salt_samples(), lw=2., alpha=1.)
Out[102]:
However, when the redshifts are not known, this is what the scatter plots look like:
In [115]:
sns.pairplot(ans.ResChar.fromSNCosmoRes(samples_z_nomodelcov).salt_samples(), diag_kind='kde', kind='scatter')
Out[115]:
Leading to posterior distributions like:
In [114]:
corner.corner(ans.ResChar.fromSNCosmoRes(samples_z_nomodelcov).salt_samples(), lw=2., alpha=1.)
Out[114]:
In [ ]: