# This code is a part of X-ray: Generate and Analyse (XGA), a module designed for the XMM Cluster Survey (XCS).
# Last modified by David J Turner (turne540@msu.edu) 20/02/2023, 14:04. Copyright (c) The Contributors
from typing import List
import numpy as np
# This set of functions is to support the MCMC fitting found in the BaseProfile class(es) - its here because
# local functions can't be pickled
[docs]def log_likelihood(theta: np.ndarray, r: np.ndarray, y: np.ndarray, y_err: np.ndarray, m_func) -> np.ndarray:
"""
Uses a simple Gaussian likelihood function, returns the logged value.
:param np.ndarray theta: The knowledge we have (think theta in Bayesian parlance) - gets fed
into the model we've chosen.
:param np.ndarray r: The radii at which we have measured profile values.
:param np.ndarray y: The values we have measured for the profile.
:param np.ndarray y_err: The uncertainties on the measured profile values.
:param m_func: The model function that is being fit to.
:return: The log-likelihood value.
:rtype: np.ndarray
"""
# Just in case something goes wrong in the model function
try:
lik = -np.sum(np.log(y_err*np.sqrt(2*np.pi)) + (((y - m_func(r, *theta))**2) / (2*y_err**2)))
except ZeroDivisionError:
lik = np.NaN
return lik
[docs]def log_prob(theta: np.ndarray, r: np.ndarray, y: np.ndarray, y_err: np.ndarray,
m_func, pr) -> np.ndarray:
"""
The combination of the log prior and log likelihood.
:param np.ndarray theta: The knowledge we have (think theta in Bayesian parlance) - gets fed
into the model we've chosen.
:param np.ndarray r: The radii at which we have measured profile values.
:param np.ndarray y: The values we have measured for the profile.
:param np.ndarray y_err: The uncertainties on the measured profile values.
:param m_func: The model function that is being fit to.
:param List pr: A list of upper and lower limits for the parameters in theta, the limits of the
uniform, uninformative priors.
:return: The log probability value.
:rtype: np.ndarray
"""
lp = log_uniform_prior(theta, pr)
if not np.isfinite(lp):
ret_val = -np.inf
else:
ret_val = lp + log_likelihood(theta, r, y, y_err, m_func)
if np.isnan(ret_val):
ret_val = -np.inf
return ret_val