In [1]:
from __future__ import print_function
import numpy as np
import sncosmo
import matplotlib as mpl
%matplotlib inline
mpl.rc('savefig', dpi=110.0)
In [2]:
data = sncosmo.read_lc('data/lc-SDSS19230.list', format='salt2')
In [3]:
# data is an astropy Table
print(data)
In [4]:
# Table metadata
print("z =", data.meta['Z_HELIO'])
print("mwebv =", data.meta['MWEBV'])
In [5]:
# Quick visualization
sncosmo.plot_lc(data);
In [6]:
model1 = sncosmo.Model(source='hsiao',
effects=[sncosmo.F99Dust()],
effect_names=['mw'],
effect_frames=['obs'])
In [7]:
print(model1)
In [8]:
# fix some values
model1.set(mwebv=data.meta['MWEBV'], z=data.meta['Z_HELIO'])
In [9]:
# fit the model to the data
result, model1 = sncosmo.fit_lc(data, model1, ['amplitude', 't0'])
In [10]:
print(result)
In [11]:
# access result components:
result.chisq, result.ndof
Out[11]:
In [12]:
sncosmo.plot_lc(data, model1, errors=result.errors);
In [13]:
# here is one row of the data table:
i = 50
print(data['Filter'][i], data['Date'][i], data['ZP'][i], data['MagSys'][i], '-->',
data['Flux'][i], '+/-', data['Fluxerr'][i])
In [14]:
# the model can predict the observed flux
model1.bandflux(data['Filter'][i], data['Date'][i], zp=data['ZP'][i], zpsys=data['MagSys'][i])
Out[14]:
In [15]:
# do it for all rows simultaneously:
predicted_flux = model1.bandflux(data['Filter'], data['Date'], zp=data['ZP'], zpsys=data['MagSys'])
predicted_flux
Out[15]:
sncosmo.fit_lc is using this function to define a likelihood or chisq, then passing that to an optimizer.
In [16]:
model2 = sncosmo.Model(source='hsiao',
effects=[sncosmo.F99Dust(), sncosmo.F99Dust()],
effect_names=['mw', 'host'],
effect_frames=['obs', 'rest'])
In [17]:
# fix known values
model2.set(mwebv=data.meta['MWEBV'], z=data.meta['Z_HELIO'])
In [18]:
print(model2)
In [19]:
# fit the model to the data
result2, model2 = sncosmo.fit_lc(data, model2, ['amplitude', 'hostebv', 't0'])
In [20]:
print(result2)
To plot the data along with the best-fit model, set the model parameters to the best-fit values and plot.
In [21]:
sncosmo.plot_lc(data, model=model2, errors=result2.errors);
In [22]:
model3 = sncosmo.Model(source='salt2',
effects=[sncosmo.F99Dust()],
effect_names=['mw'],
effect_frames=['obs'])
In [23]:
# fix known values
model3.set(mwebv=data.meta['MWEBV'], z=data.meta['Z_HELIO'])
In [24]:
# fit the model to the data
result3, model3 = sncosmo.fit_lc(data, model3, ['x0', 'x1', 'c', 't0'],
bounds={'x1': (-3., 3.), 'c':(-0.3, 0.3)})
In [25]:
print(result3)
In [26]:
sncosmo.plot_lc(data, model=model3, errors=result3.errors);
In [27]:
print(model1.get('amplitude'), " versus ", model3.get('x0'))
In [28]:
model1.source_peakmag('bessellb', 'vega')
Out[28]:
In [29]:
model3.source_peakmag('bessellb', 'vega')
Out[29]:
In [30]:
# fit the model to the data
result4, model3 = sncosmo.fit_lc(data, model3, ['x0', 'x1', 'c', 't0', 'z'],
bounds={'x1': (-3., 3.), 'c':(-0.3, 0.3), 'z':(0.25, 0.45)})
In [31]:
result4.ncall
Out[31]:
In [32]:
sncosmo.plot_lc(data, model=model3, errors=result4.errors);
In [33]:
# takes about a minute
result, model3 = sncosmo.mcmc_lc(data, model3, ['x0', 'x1', 'c', 't0', 'z'],
bounds={'x1': (-3., 3.), 'c':(-0.3, 0.3), 'z':(0.05, 0.45)},
nwalkers=50, nburn=100, nsamples=400)
In [34]:
print(result)
In [35]:
import corner
In [36]:
corner.corner(result.samples, bins=20, labels=result.vparam_names, show_titles=True);
In [ ]: