# -*- coding: UTF-8 -*-
"""
Classes defining the likelihood and the posterior probability of the model
given the data
"""
from __future__ import absolute_import
from __future__ import print_function
from __future__ import unicode_literals
import numpy as np
from celerite.modeling import Model
from scipy.stats import norm, halfcauchy
from . import io
from .passband import get_model_synmags
__all__=['WDmodel_Likelihood', 'WDmodel_Posterior', 'setup_likelihood']
[docs]def setup_likelihood(params):
"""
Setup the form of the likelihood of the data given the model.
Parameters
----------
params : dict
A parameter dictionary used to configure the
:py:class:`WDmodel_Likelihood` instance. The format of the dict is
defined by :py:func:`WDmodel.io.read_params`.
Returns
-------
lnlike : :py:class:`WDmodel.likelihood.WDmodel_Likelihood`
An instance of the likelihood function class.
"""
# parse the params to create a dictionary and init the
# WDmodel.likelihood.WDmodel_Likelihood instance
setup_args = {}
bounds = []
fixed = {}
for param in io._PARAMETER_NAMES:
setup_args[param] = params[param]['value']
bounds.append(params[param]['bounds'])
fixed[param] = params[param]['fixed']
setup_args['bounds'] = bounds
# configure the likelihood function
lnlike = WDmodel_Likelihood(**setup_args)
# freeze any parameters that we want fixed
for param, val in fixed.items():
if val:
message = "Freezing {}".format(param)
print(message)
lnlike.freeze_parameter(param)
return lnlike
class WDmodel_Likelihood(Model):
"""
Classes defining the likelihood of the model given the data
Notes
-----
Subclasses :py:class:`celerite.modeling.Model` to implement a
likelihood model for the WD atmosphere This allows parameters to be
frozen/thawed dynamically based on cmdline args/a param file.
"""
# defines the parameter names of the model
parameter_names = io._PARAMETER_NAMES
def get_value(self, spec, phot, model, covmodel, pbs, pixel_scale=1., phot_dispersion=0.):
"""
Returns the log likelihood of the model
Parameters
----------
spec : :py:class:`numpy.recarray`
The spectrum with ``dtype=[('wave', '<f8'), ('flux', '<f8'), ('flux_err', '<f8')]``
phot : None or :py:class:`numpy.recarray`
The photometry with ``dtype=[('pb', 'str'), ('mag', '<f8'), ('mag_err', '<f8')]``
model : :py:class:`WDmodel.WDmodel.WDmodel` instance
The DA White Dwarf SED model generator
covmodel : :py:class:`WDmodel.covariance.WDmodel_CovModel` instance
The parametrized model for the covariance of the spectrum ``spec``
pbs : dict
Passband dictionary containing the passbands corresponding to
``phot.pb`` and generated by :py:func:`WDmodel.passband.get_pbmodel`.
pixel_scale : float, optional
Jacobian of the transformation between wavelength in Angstrom and
pixels. In principle, this should be a vector, but virtually all
spectral reduction packages resample the spectrum onto a uniform
wavelength scale that is close to the native pixel scale of the
spectrograph. Default is ``1.``
phot_dispersion : float, optional
Excess photometric dispersion to add in quadrature with the
photometric uncertainties ``phot.mag_err``. Use if the errors are
grossly underestimated. Default is ``0.``
Returns
-------
lnlike : float
The likelihood of the model parameters :py:attr:`parameter_names`
given the data - the spectrum ``spec`` and photometry ``phot``.
"""
if phot is None:
phot_chi = 0.
mod = model._get_obs_model(self.teff, self.logg, self.av, self.fwhm,\
spec.wave, rv=self.rv, pixel_scale=pixel_scale)
else:
mod, full = model._get_full_obs_model(self.teff, self.logg, self.av, self.fwhm,\
spec.wave, rv=self.rv, pixel_scale=pixel_scale)
mod_mags = get_model_synmags(full, pbs, mu=self.mu)
phot_res = phot.mag - mod_mags.mag
phot_chi = np.sum(phot_res**2./((phot.mag_err**2.)+(phot_dispersion**2.)))
mod *= (1./(4.*np.pi*(self.dl)**2.))
res = spec.flux - mod
return covmodel.lnlikelihood(spec.wave, res, spec.flux_err, self.fsig, self.tau, self.fw) - (phot_chi/2.)
[docs]class WDmodel_Posterior(object):
"""
Classes defining the posterior probability of the model given the data
An instance of this class is used to store the data and model, and evaluate
the likelihood and prior to compute the posterior.
Parameters
----------
spec : :py:class:`numpy.recarray`
The spectrum with ``dtype=[('wave', '<f8'), ('flux', '<f8'), ('flux_err', '<f8')]``
phot : None or :py:class:`numpy.recarray`
The photometry with ``dtype=[('pb', 'str'), ('mag', '<f8'), ('mag_err', '<f8')]``
model : :py:class:`WDmodel.WDmodel.WDmodel` instance
The DA White Dwarf SED model generator
covmodel : :py:class:`WDmodel.covariance.WDmodel_CovModel` instance
The parametrized model for the covariance of the spectrum ``spec``
pbs : dict
Passband dictionary containing the passbands corresponding to
``phot.pb`` and generated by :py:func:`WDmodel.passband.get_pbmodel`.
lnlike : :py:class:`WDmodel_Likelihood` instance
Instance of the likelihood function class, such as that produced by
:py:meth:`WDmodel.likelihood.setup_likelihood`
pixel_scale : float, optional
Jacobian of the transformation between wavelength in Angstrom and
pixels. In principle, this should be a vector, but virtually all
spectral reduction packages resample the spectrum onto a uniform
wavelength scale that is close to the native pixel scale of the
spectrograph. Default is ``1.``
phot_dispersion : float, optional
Excess photometric dispersion to add in quadrature with the
photometric uncertainties ``phot.mag_err``. Use if the errors are
grossly underestimated. Default is ``0.``
Attributes
----------
spec : :py:class:`numpy.recarray`
The spectrum with ``dtype=[('wave', '<f8'), ('flux', '<f8'), ('flux_err', '<f8')]``
wave_scale : float
length of the wavelength array ``wave`` in Angstroms
phot : None or :py:class:`numpy.recarray`
The photometry with ``dtype=[('pb', 'str'), ('mag', '<f8'), ('mag_err', '<f8')]``
model : :py:class:`WDmodel.WDmodel.WDmodel` instance
The DA White Dwarf SED model generator
covmodel : :py:class:`WDmodel.covariance.WDmodel_CovModel` instance
The parametrized model for the covariance of the spectrum ``spec``
pbs : dict
Passband dictionary containing the passbands corresponding to
``phot.pb`` and generated by :py:func:`WDmodel.passband.get_pbmodel`.
_lnlike : :py:class:`WDmodel_Likelihood` instance
Instance of the likelihood function class, such as that produced by
:py:meth:`WDmodel.likelihood.setup_likelihood`
pixel_scale : float
Jacobian of the transformation between wavelength in Angstrom and
pixels. In principle, this should be a vector, but virtually all
spectral reduction packages resample the spectrum onto a uniform
wavelength scale that is close to the native pixel scale of the
spectrograph. Default is ``1.``
phot_dispersion : float, optional
Excess photometric dispersion to add in quadrature with the
photometric uncertainties ``phot.mag_err``. Use if the errors are
grossly underestimated. Default is ``0.``
p0 : dict
initial values of all the model parameters, including fixed parameters
Returns
-------
lnpost : :py:class:`WDmodel_Posterior` instance
It is this instance that is passed to the samplers/optimizers in the
:py:mod:`WDmodel.fit` module. Those methods evaluate the posterior
probability of the model parameters given the data.
Notes
-----
Wraps :py:meth:`celerite.modeling.Model.log_prior` which imposes a
boundscheck and returns ``-inf``. This is not an issue as
the samplers used in the methods in :py:mod:`WDmodel.fit`.
"""
[docs] def __init__(self, spec, phot, model, covmodel, pbs, lnlike, pixel_scale=1., phot_dispersion=0.):
self.spec = spec
self.wavescale = spec.wave.ptp()
self.phot = phot
self.model = model
self.covmodel = covmodel
self.pbs = pbs
self._lnlike = lnlike
self.pixscale = pixel_scale
self.phot_dispersion = phot_dispersion
init_p0 = lnlike.get_parameter_dict(include_frozen=True)
self.p0 = init_p0
[docs] def __call__(self, theta, prior=False, likelihood=False):
"""
Evalulates the log posterior of the model parameters given the data
Parameters
----------
theta : array-like
Vector of the non-frozen model parameters. The order of the
parameters is defined by
:py:attr:`WDmodel_Likelihood.parameter_names`.
prior : bool, optional
Only return the value of the log prior given the model parameters
likelihood : bool, optional
Only return the value of the log likelihood given the model
parameters if the prior is finite
Returns
-------
lnpost : float
the log posterior of the model parameters given the data
"""
self._lnlike.set_parameter_vector(theta)
out = self._lnprior()
if not np.isfinite(out):
return -np.inf
if prior:
return out
loglike = self._lnlike.get_value(self.spec, self.phot, self.model, self.covmodel, self.pbs,\
pixel_scale=self.pixscale, phot_dispersion=self.phot_dispersion)
if likelihood:
return loglike
out += loglike
return out
[docs] def lnlike(self, theta):
"""
Evalulates the log likelihood of the model parameters given the data.
Convenience function that can return the value of the likelihood even
if the prior is not finite unlike
:py:meth:`WDmodel_Posterior.__call__` for debugging.
Parameters
----------
theta : array-like
Vector of the non-frozen model parameters. The order of the
parameters is defined by
:py:attr:`WDmodel_Likelihood.parameter_names`.
Returns
-------
lnlike : float
the log likelihood of the model parameters given the data
"""
self._lnlike.set_parameter_vector(theta)
out = self._lnlike.get_value(self.spec, self.phot, self.model, self.covmodel, self.pbs,\
pixel_scale=self.pixscale, phot_dispersion=self.phot_dispersion)
return out
[docs] def _lnprior(self):
"""
Evalulates the log likelihood of the model parameters given the data.
Implements an lnprior function which imposes weakly informative priors
on the model parameters.
Parameters
----------
theta : array-like
Vector of the non-frozen model parameters. The order of the
parameters is defined by
:py:attr:`WDmodel_Likelihood.parameter_names`.
Returns
-------
lnprior : float
the log likelihood of the model parameters given the data
Notes
-----
The prior on ``av`` is the 'glos' prior
The prior on ``rv`` is a Gaussian with mean 3.1 and standard
deviation 0.18. This is adopted from Schlafly et al., 2014 PS1
analysis. Note that they report 3.31, but they aren't really
measuring E(B-V) with PS1. Their sigma should be consistent despite
the different filter set.
The prior on ``fsig`` and ``fw`` - the fractional amplitudes of the
non-trivial stationary and white components of the kernel used to
parametrize the covariance is half-Cauchy since we don't want it to
be less than zero
There is no explicit prior on ``tau`` i.e. a tophat prior, defined
by the bounds
The ``fwhm`` has a lower bound set at the value below which the
spectrum isn't being convolved anymore. We never run into this
bound since real spectra have physical instrumental broadening.
This prevents ``fwhm`` from going to zero for fitting poorly
simulated spectra generated from simply resampling the model grid.
The prior on all other parameters are broad Gaussians
Wraps :py:meth:`celerite.modeling.Model.log_prior` which imposes a
boundscheck and returns ``-inf``. This is not an issue as the
samplers used in the methods in :py:mod:`WDmodel.fit`.
"""
lp = self._lnlike.log_prior()
# check if out of bounds
if not np.isfinite(lp):
return -np.inf
else:
out = 0.
# put some weak priors on intrinsic WD parameters
# normal on teff
teff = self._lnlike.get_parameter('teff')
teff0 = self.p0['teff']
out += norm.logpdf(teff, teff0, 10000.)
# normal on logg
logg = self._lnlike.get_parameter('logg')
logg0 = self.p0['logg']
out += norm.logpdf(logg, logg0, 1.)
# this implements the glos prior on Av
av = self._lnlike.get_parameter('av')
avtau = 0.4
avsdelt = 0.1
wtexp = 1.
wtdelt = 0.5
sqrt2pi = np.sqrt(2.*np.pi)
normpeak = (wtdelt/sqrt2pi)/(2.*avsdelt) + wtexp/avtau
pav = (wtdelt/sqrt2pi)*np.exp((-av**2.)/(2.*avsdelt**2.))/(2.*avsdelt) +\
wtexp*np.exp(-av/avtau)/avtau
pav /= normpeak
out += np.log(pav)
# this is from the Schlafly et al analysis from PS1
rv = self._lnlike.get_parameter('rv')
out += norm.logpdf(rv, 3.1, 0.18)
# normal on dl
dl = self._lnlike.get_parameter('dl')
dl0 = self.p0['dl']
out += norm.logpdf(dl, dl0, 1000.)
fwhm = self._lnlike.get_parameter('fwhm')
# The FWHM is converted into a gaussian sigma for convolution.
# That convolution kernel is truncated at 4 standard deviations by default.
# If twice that scale is less than 1 pixel, then we're not actually modifying the data.
# This is what sets the hard lower bound on the data, not fwhm=0.
gsig = (fwhm/self.model._fwhm_to_sigma)*self.pixscale
if 8.*gsig < 1.:
return -np.inf
# normal on fwhm
fwhm0 = self.p0['fwhm']
out += norm.logpdf(fwhm, fwhm0, 8.)
# half-Cauchy on the kernel amplitudes (both stationary and white)
fsig = self._lnlike.get_parameter('fsig')
out += halfcauchy.logpdf(fsig, loc=0, scale=3)
fw = self._lnlike.get_parameter('fw')
out += halfcauchy.logpdf(fw, loc=0, scale=3)
# normal on mu
mu = self._lnlike.get_parameter('mu')
mu0 = self.p0['mu']
out += norm.logpdf(mu, mu0, 10.)
return out
[docs] def lnprior(self, theta):
"""
Evalulates the log prior of the model parameters.
Convenience function that can return the value of the prior defined to
make the interface consistent with the
:py:meth:`WDmodel_Posterior.lnlike` method. Just a thin wrapper around
:py:meth:`WDmodel_Posterior._lnprior` which is what is actually
evalulated by :py:meth:`WDmodel_Posterior.__call__`.
Parameters
----------
theta : array-like
Vector of the non-frozen model parameters. The order of the
parameters is defined by
:py:attr:`WDmodel_Likelihood.parameter_names`.
Returns
-------
lnprior : float
the log prior of the model parameters
"""
self._lnlike.set_parameter_vector(theta)
return self._lnprior()