Source code for hierarc.Likelihood.LensLikelihood.kin_likelihood

__author__ = "sibirrer"

from lenstronomy.Util import constants as const
import numpy as np


[docs] class KinLikelihood(object): """Likelihood to deal with IFU kinematics constraints with covariances in both the model and measured velocity dispersion.""" def __init__( self, z_lens, z_source, sigma_v_measurement, j_model, error_cov_measurement, error_cov_j_sqrt, normalized=True, sigma_sys_error_include=False, ): """ :param z_lens: lens redshift :param z_source: source redshift :param sigma_v_measurement: numpy array, velocity dispersion measured :param j_model: numpy array of the predicted dimensionless dispersion on the IFU's :param error_cov_measurement: covariance matrix of the measured velocity dispersions in the IFU's :param error_cov_j_sqrt: covariance matrix of sqrt(J) of the model predicted dimensionless dispersion on the IFU's :param normalized: bool, if True, returns the normalized likelihood, if False, separates the constant prefactor (in case of a Gaussian 1/(sigma sqrt(2 pi)) ) to compute the reduced chi2 statistics :param sigma_sys_error_include: bool, if True will include a systematic error in the velocity dispersion measurement (if sampled from), otherwise this sampled value is ignored. """ self._z_lens = z_lens self._j_model = np.array(j_model, dtype=float) self._sigma_v_measured = np.array(sigma_v_measurement, dtype=float) self._error_cov_measurement = np.array(error_cov_measurement, dtype=float) self._error_cov_j_sqrt = np.array(error_cov_j_sqrt, dtype=float) self.num_data = len(j_model) self._normalized = normalized self._sigma_sys_error_include = sigma_sys_error_include
[docs] def log_likelihood( self, ddt, dd, kin_scaling=None, sigma_v_sys_error=None, sigma_v_sys_offset=None, ): """ Note: kinematics + imaging data can constrain Ds/Dds. The input of Ddt, Dd is transformed here to match Ds/Dds :param ddt: time-delay distance :param dd: angular diameter distance to the deflector :param kin_scaling: array of size of the velocity dispersion measurement or None, scaling of the predicted dimensionless quantity J (proportional to sigma_v^2) of the anisotropy model in the sampling relative to the anisotropy model used to derive the prediction and covariance matrix in the init of this class. :param sigma_v_sys_error: float (optional) added error on the velocity dispersion measurement in quadrature :param sigma_v_sys_offset: float (optional) for a fractional systematic offset in the kinematic measurement such that sigma_v = sigma_v_measured * (1 + sigma_v_sys_offset) :return: log likelihood given the single lens analysis """ ds_dds = np.maximum(ddt / dd / (1 + self._z_lens), 0) if kin_scaling is None: scaling_ifu = 1 else: scaling_ifu = kin_scaling sigma_v_predict = self.sigma_v_model(ds_dds, scaling_ifu) delta = self.sigma_v_measurement_mean(sigma_v_sys_offset) - sigma_v_predict cov_error = self.cov_error_measurement( sigma_v_sys_error ) + self.cov_error_model(ds_dds, scaling_ifu) try: cov_error_inv = np.linalg.inv(cov_error) except: return -np.inf lnlikelihood = -delta.dot(cov_error_inv.dot(delta)) / 2.0 if self._normalized is True: sign_det, lndet = np.linalg.slogdet(cov_error) if sign_det < 0: raise ValueError( "error covariance matrix needs to be positive definite" ) lnlikelihood -= 1 / 2.0 * (self.num_data * np.log(2 * np.pi) + lndet) return lnlikelihood
[docs] def sigma_v_measurement_mean(self, sigma_v_sys_offset=None): """ :param sigma_v_sys_offset: float (optional) for a fractional systematic offset in the kinematic measurement such that sigma_v = sigma_v_measured * (1 + sigma_v_sys_offset) :return: corrected measured velocity dispersion """ if sigma_v_sys_offset is None: return self._sigma_v_measured else: return self._sigma_v_measured * (1 + sigma_v_sys_offset)
[docs] def sigma_v_model(self, ds_dds, kin_scaling=1): """Model predicted velocity dispersion for the IFU's. :param ds_dds: Ds/Dds :param kin_scaling: scaling of the anisotropy affecting sigma_v^2 :return: array of predicted velocity dispersions """ sigma_v_predict = np.sqrt(self._j_model * ds_dds * kin_scaling) * const.c / 1000 return sigma_v_predict
[docs] def cov_error_model(self, ds_dds, kin_scaling=1): """ :param ds_dds: Ds/Dds :param kin_scaling: scaling of the anisotropy affecting sigma_v^2 :return: covariance matrix of the error in the predicted model (from mass model uncertainties) """ scaling_matix = np.outer(np.sqrt(kin_scaling), np.sqrt(kin_scaling)) return self._error_cov_j_sqrt * scaling_matix * ds_dds * (const.c / 1000) ** 2
[docs] def cov_error_measurement(self, sigma_v_sys_error=None): """ :param sigma_v_sys_error: float (optional) added error on the velocity dispersion measurement in quadrature :return: error covariance matrix of the velocity dispersion measurements """ if self._sigma_sys_error_include and sigma_v_sys_error is not None: return self._error_cov_measurement + np.outer( self._sigma_v_measured * sigma_v_sys_error, self._sigma_v_measured * sigma_v_sys_error, ) else: return self._error_cov_measurement
[docs] def sigma_v_prediction(self, ddt, dd, kin_scaling=1): """Model prediction mean velocity dispersion vector and model prediction covariance matrix. :param ddt: time-delay distance :param dd: angular diameter distance to the deflector :param kin_scaling: array of size of the velocity dispersion measurement or None, scaling of the predicted dimensionless quantity J (proportional to sigma_v^2) of the anisotropy model in the sampling relative to the anisotropy model used to derive the prediction and covariance matrix in the init of this class. :return: model prediction mean velocity dispersion vector and model prediction covariance matrix """ ds_dds = np.maximum(ddt / dd / (1 + self._z_lens), 0) sigma_v_predict = self.sigma_v_model(ds_dds, kin_scaling) cov_error_predict = self.cov_error_model(ds_dds, kin_scaling) return sigma_v_predict, cov_error_predict
[docs] def sigma_v_measurement(self, sigma_v_sys_error=None, sigma_v_sys_offset=None): """ :param sigma_v_sys_error: float (optional) added error on the velocity dispersion measurement in quadrature :param sigma_v_sys_offset: float (optional) for a fractional systematic offset in the kinematic measurement such that sigma_v = sigma_v_measured * (1 + sigma_v_sys_offset) :return: measurement mean (vector), measurement covariance matrix """ return self.sigma_v_measurement_mean( sigma_v_sys_offset ), self.cov_error_measurement(sigma_v_sys_error)