petitRADTRANS.retrieval.log_likelihood
======================================

.. py:module:: petitRADTRANS.retrieval.log_likelihood


Classes
-------

.. autoapisummary::

   petitRADTRANS.retrieval.log_likelihood.LogLikelihood


Functions
---------

.. autoapisummary::

   petitRADTRANS.retrieval.log_likelihood.get_bic
   petitRADTRANS.retrieval.log_likelihood.get_aic
   petitRADTRANS.retrieval.log_likelihood.get_bpics


Module Contents
---------------

.. py:function:: get_bic(max_logl: float, n_data: int, n_params: int) -> float

   Calculate the Bayesian Information Criterion for a fitted model.

   References:
       Schwartz (1978), DOI 10.1214/aos/1176344136.
       Thorngren et al. (2025), DOI 10.3847/1538-4365/ae0e71.


.. py:function:: get_aic(max_logl: float, n_params: int) -> float

   Calculate the Akaike Information Criterion for a fitted model.

   References:
       Akaike (1974), DOI 10.1109/TAC.1974.1100705.
       Thorngren et al. (2025), DOI 10.3847/1538-4365/ae0e71.


.. py:function:: get_bpics(logl_samples: jax.typing.ArrayLike, n_params: int, log_weights: jax.typing.ArrayLike = None) -> float

   Calculate the simplified Bayesian Predictive Information Criterion.

   References:
       Ando (2011), DOI 10.1080/01966324.2011.10737798.
       Thorngren et al. (2025), DOI 10.3847/1538-4365/ae0e71.


.. py:class:: LogLikelihood(model_spectrum: jax.typing.ArrayLike, data_spectrum: jax.typing.ArrayLike, data_uncertainties: jax.typing.ArrayLike = None, data_covariance: jax.typing.ArrayLike = None, covariance_cholesky_factor: jax.typing.ArrayLike = None, log_likelihood_normalization: jax.typing.ArrayLike = None, n_free_parameters: int = 0)

   .. py:attribute:: model_spectrum


   .. py:attribute:: data_spectrum


   .. py:attribute:: n_data_points


   .. py:attribute:: n_free_parameters
      :value: 0



   .. py:attribute:: cholesky_factor
      :value: None



   .. py:attribute:: log_det_covariance
      :value: None



   .. py:method:: get_chi_squared() -> float


   .. py:method:: _calculate_chi_squared_covariance()


   .. py:method:: _calculate_chi_squared_uncertainties()


   .. py:method:: calculate_log_likelihood() -> float


   .. py:method:: get_reduced_chi_squared() -> float

      Calculates the reduced chi-squared value.

      Returns:
          The reduced chi-squared value. Returns jnp.nan if degrees of freedom is <= 0.



   .. py:method:: get_noise_normalization_constant() -> float

      Returns the noise-related normalization constant of the log-likelihood.
      This is -0.5 * (N*log(2pi) + log(det(Cov))) or -0.5 * sum(log(2pi*sigma_i^2)).
      This term does not depend on the model-data residuals.

      Returns:
          The pre-computed noise normalization constant.



   .. py:method:: log_likelihood(beta=None, beta_mode='multiply')


   .. py:method:: _log_likelihood(model, data, uncertainties, beta=None, beta_mode='multiply')
      :staticmethod:


      Calculate the log-likelihood between the model and the data.

      The spectrum model must be on the same wavelength grid than the data.

      From Gibson et al. 2020 (https://doi.org/10.1093/mnras/staa228). Constant terms are dropped and constant
      coefficients are set to 1.
      The 'add' beta mode comes from Line et al. 2015 (DOI 10.1088/0004-637X/807/2/183).

      Set:
          chi2(A, B) = sum(((f_i - A * m_i) / (B * sig_i)) ** 2),  implicit sum on i
      with f_i the data, m_i the model, and sig_i the uncertainty; where "i" denotes wavelength/time variation.
      Starting from (implicit product on i):
          L        = prod(1 / sqrt(2 * pi * B * sig_i ** 2) * exp(-1/2 * sum(((f_i - A * m_i) / (B * sig_i)) ** 2))),
          => L     = prod( 1 / sqrt(2 * pi * B * sig_i ** 2) * exp(-1/2 * chi2(A, B)) ),
          => ln(L) = -N/2 * ln(2 * pi) - N * ln(B) - sum(ln(sig_i)) - 1/2 * chi2(A, B).
      Dropping constant terms:
          ln(L)^* = - N * ln(B) - 1/2 * chi2(A, B).

      B can be automatically optimised by nulling the ln(L) partial derivative with respect to B.
      Using the best estimator of B instead of the true value:
          d_ln(L) / d_B = - N / B + 1 / B ** 3 * chi2(A, B=1) = 0,
          => B = sqrt( 1/N *  chi2(A, B=1)).
      Replacing:
          ln(L)^**= - N * ln(B) - 1/2 * chi2(A, B),
                  = - N * ln(sqrt(1/N * chi2(A, B=1))) - 1/2 * sum(((f_i - A * m_i) / (B * sig_i)) ** 2),
                  = - N/2 * ln(1/N * chi2(A, B=1)) - 1/2 / B ** 2 * sum(((f_i - A * m_i) / sig_i) ** 2),  B cst with i
                  = - N/2 * ln(1/N * chi2(A, B=1)) - 1/2 * N / chi2(A, B=1) * chi2(A, B=1),
                  = - N/2 * ln(1/N * chi2(A, B=1)) - N/2.
      Dropping constant terms:
          ln(L)^*** = - N/2 * ln(1/N * chi2(A, B=1)).

      Args:
          model: numpy.ndarray
              The model flux in the same units as the data.
          data: numpy.ndarray
              The data.
          uncertainties: numpy.ndarray
              The uncertainties on the data.
          beta: float, optional
              Noise scaling coefficient. If None, the noise scaling is "automatically optimised".
          beta_mode: string, optional


      Returns:
          logL : float
              The log likelihood of the model given the data.



   .. py:method:: log_likelihood2chi2(log_likelihood: float) -> float
      :staticmethod:



