|
| 1 | +import numba as nb |
| 2 | +import numpy as np |
| 3 | +from scipy.special import gammaincc, gammaincinv |
| 4 | + |
| 5 | +from preliz.distributions.distributions import Continuous |
| 6 | +from preliz.internal.distribution_helper import all_not_none, eps |
| 7 | +from preliz.internal.optimization import optimize_ml |
| 8 | +from preliz.internal.special import cdf_bounds, digamma, gammaln, ppf_bounds_cont |
| 9 | + |
| 10 | + |
| 11 | +class ScaledInverseChiSquared(Continuous): |
| 12 | + r""" |
| 13 | + Scaled Inverse Chi squared distribution. |
| 14 | +
|
| 15 | + The pdf of this distribution is |
| 16 | +
|
| 17 | + .. math:: |
| 18 | +
|
| 19 | + f(x \mid \nu, \tau^2) = |
| 20 | + \frac{(\tau^2\nu/2)^{\nu/2}}{\Gamma(\nu/2)}~ |
| 21 | + \frac{\exp\left[ \frac{-\nu \tau^2}{2 x}\right]}{x^{1+\nu/2}} |
| 22 | +
|
| 23 | +
|
| 24 | + .. plot:: |
| 25 | + :context: close-figs |
| 26 | +
|
| 27 | +
|
| 28 | + from preliz import ScaledInverseChiSquared, style |
| 29 | + style.use('preliz-doc') |
| 30 | + nus = [4., 4, 10.] |
| 31 | + tau2s = [1., 2, 1.] |
| 32 | + for nu, tau2 in zip(nus, tau2s): |
| 33 | + ScaledInverseChiSquared(nu, tau2).plot_pdf(support=(0, 5)) |
| 34 | +
|
| 35 | + ======== =============================== |
| 36 | + Support :math:`x \in [0, \infty)` |
| 37 | + Mean :math:`\nu \tau^2 / (\nu - 2)` for :math:`\nu > 2`, else :math:`\infty` |
| 38 | + Variance :math:`\frac{2 \nu^2 \tau^4}{(\nu - 2)^2 (\nu - 4)}` |
| 39 | + for :math:`\nu > 4`, else :math:`\infty` |
| 40 | + ======== =============================== |
| 41 | +
|
| 42 | +
|
| 43 | + Parameters |
| 44 | + ---------- |
| 45 | + nu : float |
| 46 | + Degrees of freedom (nu > 0). |
| 47 | + tau2 : float |
| 48 | + Scale (tau2 > 0). |
| 49 | + """ |
| 50 | + |
| 51 | + def __init__(self, nu=None, tau2=None): |
| 52 | + super().__init__() |
| 53 | + self.nu = nu |
| 54 | + self.tau2 = tau2 |
| 55 | + self.support = (0, np.inf) |
| 56 | + self._parametrization(nu, tau2) |
| 57 | + |
| 58 | + def _parametrization(self, nu=None, tau2=None): |
| 59 | + self.nu = nu |
| 60 | + self.tau2 = tau2 |
| 61 | + self.param_names = ("nu", "tau2") |
| 62 | + self.params_support = ((eps, np.inf), (eps, np.inf)) |
| 63 | + self.params = (self.nu, self.tau2) |
| 64 | + if all_not_none(nu, tau2): |
| 65 | + self._update(self.nu, self.tau2) |
| 66 | + |
| 67 | + def _update(self, nu, tau2): |
| 68 | + self.nu = np.float64(nu) |
| 69 | + self.tau2 = np.float64(tau2) |
| 70 | + self.params = (self.nu, self.tau2) |
| 71 | + self.is_frozen = True |
| 72 | + |
| 73 | + def pdf(self, x): |
| 74 | + x = np.asarray(x) |
| 75 | + return np.exp(nb_logpdf(x, self.nu, self.tau2)) |
| 76 | + |
| 77 | + def cdf(self, x): |
| 78 | + x = np.asarray(x) |
| 79 | + return nb_cdf(x, self.nu, self.tau2) |
| 80 | + |
| 81 | + def ppf(self, q): |
| 82 | + q = np.asarray(q) |
| 83 | + return nb_ppf(q, self.nu, self.tau2) |
| 84 | + |
| 85 | + def logpdf(self, x): |
| 86 | + return nb_logpdf(x, self.nu, self.tau2) |
| 87 | + |
| 88 | + def _neg_logpdf(self, x): |
| 89 | + return nb_neg_logpdf(x, self.nu, self.tau2) |
| 90 | + |
| 91 | + def entropy(self): |
| 92 | + return nb_entropy(self.nu, self.tau2) |
| 93 | + |
| 94 | + def mean(self): |
| 95 | + return np.where(self.nu > 2, self.nu * self.tau2 / (self.nu - 2), np.inf) |
| 96 | + |
| 97 | + def mode(self): |
| 98 | + return self.nu * self.tau2 / (self.nu + 2) |
| 99 | + |
| 100 | + def median(self): |
| 101 | + return self.ppf(0.5) |
| 102 | + |
| 103 | + def var(self): |
| 104 | + return np.where( |
| 105 | + self.nu > 4, |
| 106 | + 2 * self.nu**2 * self.tau2**2 / ((self.nu - 2) ** 2 * (self.nu - 4)), |
| 107 | + np.inf, |
| 108 | + ) |
| 109 | + |
| 110 | + def std(self): |
| 111 | + return self.var() ** 0.5 |
| 112 | + |
| 113 | + def skewness(self): |
| 114 | + return np.where(self.nu > 6, 4 / (self.nu - 6) * np.sqrt(2 * (self.nu - 4)), np.inf) |
| 115 | + |
| 116 | + def kurtosis(self): |
| 117 | + return np.where( |
| 118 | + self.nu > 8, (12 * (5 * self.nu - 22)) / ((self.nu - 6) * (self.nu - 8)), np.inf |
| 119 | + ) |
| 120 | + |
| 121 | + def rvs(self, size=None, random_state=None): |
| 122 | + random_state = np.random.default_rng(random_state) |
| 123 | + return (self.nu * self.tau2) / random_state.chisquare(self.nu, size) |
| 124 | + |
| 125 | + def _fit_moments(self, mean, sigma): |
| 126 | + cv2 = sigma**2 / mean**2 |
| 127 | + nu_hat = 4 + 2 / cv2 |
| 128 | + tau2_hat = mean * (nu_hat - 2) / nu_hat |
| 129 | + self._update(nu_hat, tau2_hat) |
| 130 | + |
| 131 | + def _fit_mle(self, sample): |
| 132 | + optimize_ml(self, sample) |
| 133 | + |
| 134 | + |
| 135 | +# @nb.njit(cache=True) |
| 136 | +def nb_cdf(x, nu, tau2): |
| 137 | + h_nu = nu / 2 |
| 138 | + return cdf_bounds(gammaincc(h_nu, h_nu * tau2 / x), x, 0, np.inf) |
| 139 | + |
| 140 | + |
| 141 | +# @nb.njit(cache=True) |
| 142 | +def nb_ppf(q, nu, tau2): |
| 143 | + h_nu = nu / 2 |
| 144 | + vals = h_nu * tau2 / gammaincinv(h_nu, 1 - q) |
| 145 | + return ppf_bounds_cont(vals, q, 0, np.inf) |
| 146 | + |
| 147 | + |
| 148 | +@nb.vectorize(nopython=True, cache=True) |
| 149 | +def nb_logpdf(x, nu, tau2): |
| 150 | + if x < 0: |
| 151 | + return -np.inf |
| 152 | + else: |
| 153 | + h_nu = nu / 2 |
| 154 | + return ( |
| 155 | + -(np.log(x) * (h_nu + 1)) |
| 156 | + - (h_nu * tau2) / x |
| 157 | + + np.log(tau2) * h_nu |
| 158 | + - gammaln(h_nu) |
| 159 | + + np.log(h_nu) * h_nu |
| 160 | + ) |
| 161 | + |
| 162 | + |
| 163 | +@nb.njit(cache=True) |
| 164 | +def nb_neg_logpdf(x, nu, tau2): |
| 165 | + return (-nb_logpdf(x, nu, tau2)).sum() |
| 166 | + |
| 167 | + |
| 168 | +@nb.njit(cache=True) |
| 169 | +def nb_entropy(nu, tau2): |
| 170 | + h_nu = nu / 2 |
| 171 | + return h_nu + np.log(h_nu * tau2) + gammaln(h_nu) - (1 + h_nu) * digamma(h_nu) |
0 commit comments