Source code for rework_pysatl_mpest.distributions.cauchy

"""Module providing Cauchy distribution class"""

__author__ = "Maksim Pastukhov"
__copyright__ = "Copyright (c) 2025 PySATL project"
__license__ = "SPDX-License-Identifier: MIT"

import numpy as np
from numpy import float64
from scipy.stats import cauchy

from ..core import Parameter
from .continuous_dist import ContinuousDistribution


[docs] class Cauchy(ContinuousDistribution): """Class for the two-parameter cauchy distribution. Parameters ---------- loc : float Location parameter. Can be any real number. scale : float Scale parameter (gamma). Must be positive. Attributes ---------- loc : float Location parameter. scale : float Scale parameter. Methods ------- .. autosummary:: :toctree: generated/ ppf pdf lpdf log_gradients generate """ PARAM_LOC = "loc" PARAM_SCALE = "scale" loc = Parameter() scale = Parameter(lambda x: x > 0.0, "Scale parameter should be positive") def __init__(self, loc: float, scale: float): super().__init__() self.loc = loc self.scale = scale @property def name(self) -> str: return "Cauchy" @property def params(self) -> set[str]: return {self.PARAM_LOC, self.PARAM_SCALE}
[docs] def pdf(self, X): """Probability density function (PDF). The PDF for the two-parameter cauchy distribution is: .. math:: f(x | \\alpha, \\beta) = \\frac{1}{\\pi \\cdot \\beta \\cdot (1 +(\\frac{(x - \\alpha)}{\\beta})^2)} where :math:`\\beta` is the scale parameter and :math:`\\alpha` is the location parameter. Parameters ---------- X : ArrayLike The input data points at which to evaluate the PDF. Returns ------- NDArray[np.float64] The PDF values corresponding to each point in :attr:`X`. """ X = np.asarray(X, dtype=float64) return 1.0 / (np.pi * self.scale * (1.0 + ((X - self.loc) / self.scale) ** 2))
[docs] def ppf(self, P): """Percent Point Function (PPF) or quantile function. The PPF for the two-parameter cauchy distribution is: .. math:: Q(p | \\alpha, \\beta) = \\alpha + \\beta \\cdot \\tan(\\pi \\cdot (p - 0.5)) where :math:`\\alpha` is the location parameter and :math:`\\beta` is the scale parameter. Parameters ---------- P : ArrayLike The probability values (between 0 and 1) at which to evaluate the PPF. Returns ------- NDArray[np.float64] The PPF values corresponding to each probability in :attr:`P`. """ P = np.asarray(P, dtype=float64) return np.where( (P >= 0) & (P <= 1), np.where( (P == 0) | (P == 1), np.where(P == 1, np.inf, -np.inf), self.loc + self.scale * np.tan(np.pi * (P - 0.5)), ), np.nan, )
[docs] def lpdf(self, X): """Log of the Probability Density Function (LPDF). The log-PDF for the two-parameter cauchy distribution is: .. math:: \\ln f(x | \\alpha, \\beta) = \\ln(\\gamma) - \\ln(\\pi \\cdot ((x - \\alpha)^2 + \\beta^2)) where :math:`\\alpha` is the location parameter and :math:`\\beta` is the scale parameter. Parameters ---------- X : ArrayLike The input data points at which to evaluate the LPDF. Returns ------- NDArray[np.float64] The log-PDF values corresponding to each point in :attr:`X`. """ X = np.asarray(X, dtype=float64) return np.log(1.0) - np.log(np.pi * self.scale * (1.0 + ((X - self.loc) / self.scale) ** 2))
def _dlog_loc(self, X): """Partial derivative of the lpdf w.r.t. the :attr:`loc` parameter. The derivative is non-zero only for `X >= loc`. .. math:: \\frac{\\partial \\ln f(x | \\alpha, \\beta)}{\\partial \\alpha} = \\frac{2 \\cdot x - 2 \\cdot \\alpha}{\\beta^2 + x^2 - 2 \\cdot \\alpha \\cdot x + \\alpha^2} where :math:`\\alpha` is the location parameter and :math:`\\beta` is the scale parameter. Parameters ---------- X : ArrayLike The input data points. Returns ------- NDArray[np.float64] The gradient of the lpdf with respect to :attr:`loc` for each point in ::attr`X`. """ X = np.asarray(X, dtype=float64) return np.where( self.loc <= X, (2 * X - 2 * self.loc) / (self.scale**2 + X**2 - 2 * self.loc * X + self.loc**2), 0.0 ) def _dlog_scale(self, X): """Partial derivative of the lpdf w.r.t. the :attr:`scale` parameter. The derivative is non-zero only for `X >= loc`. .. math:: \\frac{\\partial \\ln f(x | \\alpha, \\beta)}{\\partial \\alpha} = \\frac{-\\beta^2 + x^2 - 2 \\cdot \\alpha \\cdot x + \\alpha^2}{\\beta^3 + \\beta \\cdot (x^2) - 2 \\cdot \\alpha \\cdot \\beta \\cdot x + \\beta \\cdot \\alpha^2} where :math:`\\alpha` is the location and :math:`\\beta` is the scale. Parameters ---------- X : ArrayLike The input data points. Returns ------- NDArray[np.float64] The gradient of the lpdf with respect to :attr:`rate` for each point in :attr:`X`. """ X = np.asarray(X, dtype=float64) return np.where( self.loc <= X, (-(self.scale**2) + X**2 - 2 * self.loc * X + self.loc**2) / (self.scale**3 + self.scale * (X**2) - 2 * self.loc * self.scale * X + self.scale * self.loc**2), 0.0, )
[docs] def log_gradients(self, X): """Calculates the gradients of the log-PDF w.r.t. its parameters. The gradients are computed for the parameters that are not fixed. Parameters ---------- X : ArrayLike The input data points at which to calculate the gradients. Returns ------- NDArray[np.float64] An array where each row corresponds to a data point in :attr:`X` and each column corresponds to the gradient with respect to a specific optimizable parameter. The order of columns corresponds to the sorted order of :attr:`self.params_to_optimize`. """ X = np.asarray(X, dtype=float64) gradient_calculators = { self.PARAM_LOC: self._dlog_loc, self.PARAM_SCALE: self._dlog_scale, } optimizable_params = sorted(list(self.params_to_optimize)) if not optimizable_params: return np.empty((len(X), 0)) gradients = [gradient_calculators[param](X) for param in optimizable_params] return np.stack(gradients, axis=1)
[docs] def generate(self, size: int): """Generates random samples from the distribution. Parameters ---------- size : int The number of random samples to generate. Returns ------- NDArray[np.float64] A NumPy array containing the generated samples. """ return np.asarray(cauchy.rvs(loc=self.loc, scale=self.scale, size=size), dtype=float64)
def __repr__(self) -> str: """Returns a string representation of the object. Returns ------- str A string that can be used to recreate the object, e.g., "Cauchy(loc=0.0, scale=2.0)". """ return f"{self.__class__.__name__}(loc={self.loc}, scale={self.scale})"