rv/dist/
inv_chi_squared.rs

1//! Χ<sup>-2</sup> over x in (0, ∞)
2#[cfg(feature = "serde1")]
3use serde::{Deserialize, Serialize};
4
5use crate::impl_display;
6use crate::misc::ln_gammafn;
7use crate::traits::{
8    Cdf, ContinuousDistr, HasDensity, Kurtosis, Mean, Mode, Parameterized,
9    Sampleable, Scalable, Shiftable, Skewness, Support, Variance,
10};
11use rand::Rng;
12use special::Gamma;
13use std::f64::consts::LN_2;
14use std::fmt;
15use std::sync::OnceLock;
16
17/// [Χ<sup>-2</sup> distribution](https://en.wikipedia.org/wiki/Inverse-chi-squared_distribution)
18/// Χ<sup>-2</sup>(v).
19///
20/// # Example
21///
22/// ```
23/// use rv::prelude::*;
24///
25/// let ix2 = InvChiSquared::new(2.0).unwrap();
26/// ```
27///
28/// Parameters for the Inverse Chi-squared distribution
29#[derive(Debug, Clone, PartialEq)]
30#[cfg_attr(feature = "serde1", derive(Serialize, Deserialize))]
31#[cfg_attr(feature = "serde1", serde(rename_all = "snake_case"))]
32pub struct InvChiSquaredParameters {
33    /// Degrees of freedom in (0, ∞)
34    pub v: f64,
35}
36
37#[derive(Debug, Clone)]
38#[cfg_attr(feature = "serde1", derive(Serialize, Deserialize))]
39#[cfg_attr(feature = "serde1", serde(rename_all = "snake_case"))]
40pub struct InvChiSquared {
41    /// Degrees of freedom in (0, ∞)
42    v: f64,
43    // ln( 2^{-v/2} / gamma(v/2))
44    #[cfg_attr(feature = "serde1", serde(skip))]
45    ln_f_const: OnceLock<f64>,
46}
47
48impl Parameterized for InvChiSquared {
49    type Parameters = InvChiSquaredParameters;
50
51    fn emit_params(&self) -> Self::Parameters {
52        Self::Parameters { v: self.v() }
53    }
54
55    fn from_params(params: Self::Parameters) -> Self {
56        Self::new_unchecked(params.v)
57    }
58}
59
60impl PartialEq for InvChiSquared {
61    fn eq(&self, other: &InvChiSquared) -> bool {
62        self.v == other.v
63    }
64}
65
66crate::impl_shiftable!(InvChiSquared);
67crate::impl_scalable!(InvChiSquared);
68
69#[derive(Debug, Clone, PartialEq)]
70#[cfg_attr(feature = "serde1", derive(Serialize, Deserialize))]
71#[cfg_attr(feature = "serde1", serde(rename_all = "snake_case"))]
72pub enum InvChiSquaredError {
73    /// v parameter is less than or equal to zero
74    VTooLow { v: f64 },
75    /// v parameter is infinite or NaN
76    VNotFinite { v: f64 },
77}
78
79impl InvChiSquared {
80    /// Create a new Inverse Chi-squared distribution
81    ///
82    /// # Arguments
83    /// - v: Degrees of freedom in (0, ∞)
84    #[inline]
85    pub fn new(v: f64) -> Result<Self, InvChiSquaredError> {
86        if v <= 0.0 {
87            Err(InvChiSquaredError::VTooLow { v })
88        } else if !v.is_finite() {
89            Err(InvChiSquaredError::VNotFinite { v })
90        } else {
91            Ok(InvChiSquared {
92                v,
93                ln_f_const: OnceLock::new(),
94            })
95        }
96    }
97
98    /// Create a new `InvChiSquared` without checking whether the parameters are
99    /// valid.
100    #[inline(always)]
101    #[must_use]
102    pub fn new_unchecked(v: f64) -> Self {
103        InvChiSquared {
104            v,
105            ln_f_const: OnceLock::new(),
106        }
107    }
108
109    /// Get the degrees of freedom, `v`.
110    ///
111    /// # Example
112    ///
113    /// ```rust
114    /// # use rv::dist::InvChiSquared;
115    /// let ix2 = InvChiSquared::new(1.2).unwrap();
116    /// assert_eq!(ix2.v(), 1.2);
117    /// ```
118    #[inline(always)]
119    pub fn v(&self) -> f64 {
120        self.v
121    }
122
123    /// Set the degrees of freedom, `k`.
124    ///
125    /// # Example
126    ///
127    /// ```rust
128    /// # use rv::dist::InvChiSquared;
129    /// let mut ix2 = InvChiSquared::new(1.2).unwrap();
130    ///
131    /// ix2.set_v(2.2).unwrap();
132    /// assert_eq!(ix2.v(), 2.2);
133    /// ```
134    ///
135    /// Will error given invalid values.
136    ///
137    /// ```rust
138    /// # use rv::dist::InvChiSquared;
139    /// # let mut ix2 = InvChiSquared::new(1.2).unwrap();
140    /// assert!(ix2.set_v(2.2).is_ok());
141    /// assert!(ix2.set_v(0.0).is_err());
142    /// assert!(ix2.set_v(-1.0).is_err());
143    /// assert!(ix2.set_v(f64::NAN).is_err());
144    /// assert!(ix2.set_v(f64::INFINITY).is_err());
145    /// ```
146    #[inline]
147    pub fn set_v(&mut self, v: f64) -> Result<(), InvChiSquaredError> {
148        if !v.is_finite() {
149            Err(InvChiSquaredError::VNotFinite { v })
150        } else if v > 0.0 {
151            self.set_v_unchecked(v);
152            Ok(())
153        } else {
154            Err(InvChiSquaredError::VTooLow { v })
155        }
156    }
157
158    #[inline(always)]
159    pub fn set_v_unchecked(&mut self, v: f64) {
160        self.v = v;
161        self.ln_f_const = OnceLock::new();
162    }
163
164    /// Get ln( 2^{-v/2} / Gamma(v/2) )
165    #[inline]
166    fn ln_f_const(&self) -> f64 {
167        *self.ln_f_const.get_or_init(|| {
168            let v2 = self.v / 2.0;
169            (-v2).mul_add(LN_2, -ln_gammafn(v2))
170        })
171    }
172}
173
174impl From<&InvChiSquared> for String {
175    fn from(ix2: &InvChiSquared) -> String {
176        format!("χ⁻²({})", ix2.v)
177    }
178}
179
180impl_display!(InvChiSquared);
181
182macro_rules! impl_traits {
183    ($kind:ty) => {
184        impl HasDensity<$kind> for InvChiSquared {
185            fn ln_f(&self, x: &$kind) -> f64 {
186                let x64 = f64::from(*x);
187                let z = self.ln_f_const();
188                (-self.v / 2.0 - 1.0).mul_add(x64.ln(), z) - (2.0 * x64).recip()
189            }
190        }
191
192        impl Sampleable<$kind> for InvChiSquared {
193            fn draw<R: Rng>(&self, rng: &mut R) -> $kind {
194                let x2 = rand_distr::ChiSquared::new(self.v).unwrap();
195                let x_inv: f64 = rng.sample(x2);
196                x_inv.recip() as $kind
197            }
198        }
199
200        impl Support<$kind> for InvChiSquared {
201            fn supports(&self, x: &$kind) -> bool {
202                *x > 0.0 && x.is_finite()
203            }
204        }
205
206        impl ContinuousDistr<$kind> for InvChiSquared {}
207
208        impl Mean<$kind> for InvChiSquared {
209            fn mean(&self) -> Option<$kind> {
210                if self.v > 2.0 {
211                    Some(1.0 / (self.v as $kind - 2.0))
212                } else {
213                    None
214                }
215            }
216        }
217
218        impl Mode<$kind> for InvChiSquared {
219            fn mode(&self) -> Option<$kind> {
220                Some((1.0 / (self.v + 2.0)) as $kind)
221            }
222        }
223
224        impl Variance<$kind> for InvChiSquared {
225            fn variance(&self) -> Option<$kind> {
226                if self.v > 4.0 {
227                    let denom =
228                        (self.v - 2.0) * (self.v - 2.0) * (self.v - 4.0);
229                    Some((2.0 / denom) as $kind)
230                } else {
231                    None
232                }
233            }
234        }
235
236        impl Cdf<$kind> for InvChiSquared {
237            fn cdf(&self, x: &$kind) -> f64 {
238                let x64 = f64::from(*x);
239                if x64 <= 0.0 {
240                    0.0
241                } else {
242                    1.0 - (2.0 * x64).recip().inc_gamma(self.v / 2.0)
243                }
244            }
245        }
246    };
247}
248
249impl Skewness for InvChiSquared {
250    fn skewness(&self) -> Option<f64> {
251        if self.v > 6.0 {
252            let v = self.v;
253            Some(4.0 / (v - 6.0) * (2.0 * (v - 4.0)).sqrt())
254        } else {
255            None
256        }
257    }
258}
259
260impl Kurtosis for InvChiSquared {
261    fn kurtosis(&self) -> Option<f64> {
262        if self.v > 8.0 {
263            let v = self.v;
264            Some(12.0 * 5.0_f64.mul_add(v, -22.0) / ((v - 6.0) * (v - 8.0)))
265        } else {
266            None
267        }
268    }
269}
270
271impl_traits!(f64);
272impl_traits!(f32);
273
274impl std::error::Error for InvChiSquaredError {}
275
276impl fmt::Display for InvChiSquaredError {
277    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
278        match self {
279            Self::VTooLow { v } => {
280                write!(f, "v ({v}) must be greater than zero")
281            }
282            Self::VNotFinite { v } => write!(f, "v ({v}) must be finite"),
283        }
284    }
285}
286
287#[cfg(test)]
288mod test {
289    use super::*;
290
291    use crate::dist::{Gamma, InvGamma};
292    use crate::misc::ks_test;
293    use crate::test_basic_impls;
294    use std::f64;
295
296    const TOL: f64 = 1E-12;
297    const KS_PVAL: f64 = 0.2;
298    const N_TRIES: usize = 5;
299
300    test_basic_impls!(f64, InvChiSquared, InvChiSquared::new(3.2).unwrap());
301
302    #[test]
303    fn new() {
304        let ix2 = InvChiSquared::new(2.3).unwrap();
305        assert::close(ix2.v, 2.3, TOL);
306    }
307
308    #[test]
309    fn new_should_reject_v_leq_zero() {
310        assert!(InvChiSquared::new(f64::MIN_POSITIVE).is_ok());
311        assert!(InvChiSquared::new(0.0).is_err());
312        assert!(InvChiSquared::new(-f64::MIN_POSITIVE).is_err());
313        assert!(InvChiSquared::new(-1.0).is_err());
314    }
315
316    #[test]
317    fn new_should_reject_non_finite_v() {
318        assert!(InvChiSquared::new(f64::INFINITY).is_err());
319        assert!(InvChiSquared::new(-f64::NAN).is_err());
320        assert!(InvChiSquared::new(f64::NEG_INFINITY).is_err());
321    }
322
323    #[test]
324    fn mean_is_defined_for_v_gt_2() {
325        {
326            let m: Option<f64> = InvChiSquared::new_unchecked(0.1).mean();
327            assert!(m.is_none());
328        }
329        {
330            let m: Option<f64> = InvChiSquared::new_unchecked(2.0).mean();
331            assert!(m.is_none());
332        }
333        {
334            let m: Option<f64> = InvChiSquared::new_unchecked(2.000_001).mean();
335            assert!(m.is_some());
336        }
337    }
338
339    #[test]
340    fn mean_values() {
341        {
342            let m: f64 = InvChiSquared::new_unchecked(2.1).mean().unwrap();
343            assert::close(m, 10.0, TOL);
344        }
345        {
346            let m: f64 = InvChiSquared::new_unchecked(4.0).mean().unwrap();
347            assert::close(m, 0.5, TOL);
348        }
349    }
350
351    #[test]
352    fn variance_is_defined_for_v_gt_4() {
353        {
354            let m: Option<f64> = InvChiSquared::new_unchecked(0.1).variance();
355            assert!(m.is_none());
356        }
357        {
358            let m: Option<f64> = InvChiSquared::new_unchecked(4.0).variance();
359            assert!(m.is_none());
360        }
361        {
362            let m: Option<f64> =
363                InvChiSquared::new_unchecked(4.000_001).variance();
364            assert!(m.is_some());
365        }
366    }
367
368    #[test]
369    fn variance() {
370        let v: f64 = InvChiSquared::new_unchecked(6.0).variance().unwrap();
371        assert::close(v, 0.0625, TOL);
372    }
373
374    #[test]
375    fn skewness() {
376        let v: f64 = InvChiSquared::new_unchecked(12.0).skewness().unwrap();
377        assert::close(v, 16.0 / 6.0, TOL);
378    }
379
380    #[test]
381    fn kurtosis() {
382        let v: f64 = InvChiSquared::new_unchecked(12.0).kurtosis().unwrap();
383        assert::close(v, 19.0, TOL);
384    }
385
386    #[test]
387    fn pdf_agrees_with_inv_gamma_special_case() {
388        let mut rng = rand::rng();
389        let v_prior = Gamma::new_unchecked(2.0, 1.0);
390
391        for v in v_prior.sample_stream(&mut rng).take(1000) {
392            let ix2 = InvChiSquared::new(v).unwrap();
393            let igam = InvGamma::new(v / 2.0, 0.5).unwrap();
394
395            for x in &[0.1_f64, 1.0_f64, 14.2_f64] {
396                assert::close(ix2.ln_f(x), igam.ln_f(x), TOL);
397            }
398        }
399    }
400
401    #[test]
402    fn cdf_limits_are_0_and_1() {
403        let ix2 = InvChiSquared::new(2.5).unwrap();
404        assert::close(ix2.cdf(&1e-16), 0.0, TOL);
405        assert::close(ix2.cdf(&1E16), 1.0, TOL);
406    }
407
408    #[test]
409    fn cdf_agrees_with_inv_gamma_special_case() {
410        let mut rng = rand::rng();
411        let v_prior = Gamma::new_unchecked(2.0, 1.0);
412
413        for v in v_prior.sample_stream(&mut rng).take(1000) {
414            let ix2 = InvChiSquared::new(v).unwrap();
415            let igam = InvGamma::new(v / 2.0, 0.5).unwrap();
416
417            for x in &[0.1_f64, 1.0_f64, 14.2_f64] {
418                assert::close(ix2.cdf(x), igam.cdf(x), TOL);
419            }
420        }
421    }
422
423    #[test]
424    fn draw_agrees_with_cdf() {
425        let mut rng = rand::rng();
426        let ix2 = InvChiSquared::new(1.2).unwrap();
427        let cdf = |x: f64| ix2.cdf(&x);
428
429        // test is flaky, try a few times
430        let passes = (0..N_TRIES).fold(0, |acc, _| {
431            let xs: Vec<f64> = ix2.sample(1000, &mut rng);
432            let (_, p) = ks_test(&xs, cdf);
433            if p > KS_PVAL { acc + 1 } else { acc }
434        });
435
436        assert!(passes > 0);
437    }
438
439    #[test]
440    fn emit_and_from_params_are_identity() {
441        let dist_a = InvChiSquared::new(1.5).unwrap();
442        let dist_b = InvChiSquared::from_params(dist_a.emit_params());
443        assert_eq!(dist_a, dist_b);
444    }
445}
rv/dist/inv_chi_squared.rs

rv/dist/
inv_chi_squared.rs
