notebook.community

Edit and run



In [1]:

    
import numpy as np



In [2]:

    
import sncosmo



In [3]:

    
import os



In [4]:

    
%matplotlib inline
import matplotlib.pyplot as plt



In [5]:

    
import analyzeSN as ans
from analyzeSN import LightCurve, ResChar



In [6]:

    
print sncosmo.__version__









    



1.5.dev722

Read in light Curves

The light curves for JLA are available from the JLA page. However a couple of them have been added to the example_data directory for easy demonstration.



In [7]:

    
sdssfname = os.path.join(ans.example_data, 'lc-SDSS14331.list')
snlsfname = os.path.join(ans.example_data, 'lc-03D1au.list')



In [8]:

    
os.path.exists(snlsfname)









    Out[8]:





True



In [9]:

    
sdsslc = LightCurve.fromSALTFormat(sdssfname)
snlslc = LightCurve.fromSALTFormat(snlsfname)

Light Curve Plots



In [10]:

    
fig_snls = sncosmo.plot_lc(snlslc.snCosmoLC())



In [11]:

    
fig_sdss = sncosmo.plot_lc(sdsslc.snCosmoLC())

Inference of Light Curve Model Parameters



In [12]:

    
import statsmodels



In [13]:

    
from scipy import stats



In [14]:

    
dust = sncosmo.CCM89Dust()
model = sncosmo.Model(source='salt2', effects=[dust], effect_names=['mw'], effect_frames=['obs'])
#model.set(mwebv=['MWEBV'])



In [22]:

    
model.set(mwebv=snlslc.props['MWEBV'])
model.set(z=snlslc.props['Redshift'])
snls_fitres = sncosmo.mcmc_lc(snlslc.snCosmoLC(), model=model,
                              bounds=dict(#z=(0.01,1.05), 
                                     c=(-0.5, 0.5), x1=(-5.,5.)),
                              vparam_names=['x0', 'x1', 'c', 't0'],
                              modelcov=True)



In [30]:

    
snls_fitres_z = sncosmo.mcmc_lc(snlslc.snCosmoLC(), model=model, bounds=dict(z=(0.05, 0.8), c=(-0.5, 0.5), x1=(-5., 5.)),
                        vparam_names=['x0', 'x1', 'c', 't0', 'z'],
                        priors=dict(z=lambda x: stats.norm.pdf(x,0.5, 0.2)), 
                        modelcov=True)



In [31]:

    
mcmc_reschar = ResChar.fromSNCosmoRes(snls_fitres)
mcmc_reschar_z = ResChar.fromSNCosmoRes(snls_fitres_z)



In [29]:

    
fig, ax = plt.subplots()
mcmc_reschar.salt_samples().mu.hist(histtype='step',
                                    bins=np.arange(10., 14.,0.05),
                                    ax=ax, normed=1)
mcmc_reschar_z.salt_samples().mu.hist(histtype='step',
                                      bins=np.arange(10., 14.,0.05),
                                      ax=ax, normed=1, color='r')


#x = np.arange(0.1, 0.4, 0.001)
#ax.plot(x, 120. *stats.norm.pdf(x, 0.213, 0.4), 'r--')









    Out[29]:





<matplotlib.axes._subplots.AxesSubplot at 0x10cd52890>



In [45]:

    
fig_z, ax_z = plt.subplots()
x = np.arange(0.1, 1.4, 0.001)
mcmc_reschar_z.salt_samples().z.hist(histtype='step', 
                                     bins=x,
                                     normed=1,
                                     ax=ax_z)
ax_z.plot(x, stats.norm.pdf(x, 0.5, 0.3), 'r-')
ax_z.axvline(snlslc.props['Redshift'], **dict(lw=2, color='k'))









    Out[45]:





<matplotlib.lines.Line2D at 0x10f076350>

Scratch



In [ ]:

    
import seaborn as sns
sns.set_style('whitegrid')



In [149]:

    
sns.pairplot(data=mcmc_reschar.salt_samples())









    Out[149]:





<seaborn.axisgrid.PairGrid at 0x10f650bd0>



In [131]:

    
ResChar.fromSNCosmoRes(fitres).parameters









    Out[131]:





z            0.851218
t0       54021.692575
x0           0.000027
x1           3.054032
c           -0.352032
mwebv        0.022100
mwr_v        3.100000
dtype: float64



In [132]:

    
fig = sncosmo.plot_lc(lco.snCosmoLC(), model=fitres[1])



In [15]:

    
fig









    Out[15]:

Setting model parameters from JLA_paramslc.txt by hand



In [16]:

    
model.set(x1=0.141351)
model.set(z=0.213)
model.set(c=0.0142)



In [17]:

    
print model









    



<Model at 0x10a48b410>
source:
  class      : SALT2Source
  name       : 'salt2'
  version    : 2.4
  phases     : [-20, .., 50] days
  wavelengths: [2000, .., 9200] Angstroms
effect (name='mw' frame='obs'):
  class           : CCM89Dust
  wavelength range: [1000, 33333.3] Angstroms
parameters:
  z     = 0.21299999999999999
  t0    = 0.0
  x0    = 1.0
  x1    = 0.141351
  c     = 0.014200000000000001
  mwebv = 0.022100000000000002
  mwr_v = 3.1000000000000001



In [20]:

    
fitres_JLAparams = sncosmo.fit_lc(lco.snCosmoLC(), model=model, vparam_names=['x0', 't0'], modelcov=True)



In [21]:

    
fitres_JLAparams[0].parameters









    Out[21]:





array([  2.13000000e-01,   5.40134824e+04,   5.79360360e-05,
         1.41351000e-01,   1.42000000e-02,   2.21000000e-02,
         3.10000000e+00])



In [22]:

    
fig = sncosmo.plot_lc(lco.snCosmoLC(), model=(fitres[1], fitres_JLAparams[1]))



In [23]:

    
fig









    Out[23]:



In [20]:

    
mcmcres_JLAparams = sncosmo.mcmc_lc(lco.snCosmoLC(), model=model, vparam_names=['x0', 't0'], modelcov=True)



In [21]:

    
ResChar.fromSNCosmoRes(mcmcres_JLAparams).parameters









    Out[21]:





z            0.213000
t0       54013.470555
x0           0.000058
x1           0.141351
c            0.014200
mwebv        0.022100
mwr_v        3.100000
dtype: float64



In [22]:

    
sncosmo.plot_lc(lco.snCosmoLC(), model=[fitres[1], fitres_JLAparams[1], mcmcres_JLAparams[1]], color='k')









    Out[22]:



In [27]:

    
fitres_JLAparams = sncosmo.fit_lc(lco.snCosmoLC(), model=model, vparam_names=['x0', 't0'], modelcov=False, minsnr=5.)



In [28]:

    
fitres_JLAparams[0].parameters









    Out[28]:





array([  2.13000000e-01,   5.40134475e+04,   5.83488019e-05,
         1.41351000e-01,   1.42000000e-02,   2.21000000e-02,
         3.10000000e+00])



In [ ]: