Source code for WDmodel.likelihood

# -*- 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()