1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126
use crate::{ consts::PI, prelude::*, }; use rand::Rng; use spaces::real::Reals; use std::fmt; pub use crate::params::Loc; new_dist!(StudentT<Loc<f64>>); macro_rules! get_nu { ($self:ident) => { ($self.0).0 } } impl Distribution for StudentT { type Support = Reals; type Params = Loc<f64>; fn support(&self) -> Reals { Reals } fn params(&self) -> Loc<f64> { self.0 } fn cdf(&self, x: &f64) -> Probability { use special_fun::FloatSpecial; let nu = get_nu!(self); let np1o2 = (nu + 1.0) / 2.0; let hyp2f1 = 0.5f64.hyp2f1(np1o2, 3.0 / 2.0, -x * x / nu); Probability::new_unchecked( 0.5 + x * np1o2.gamma() * hyp2f1 / (nu * PI).sqrt() * (nu / 2.0).gamma() ) } fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> f64 { use rand_distr::Distribution as _; rand_distr::StudentT::new(get_nu!(self)).unwrap().sample(rng) } } impl ContinuousDistribution for StudentT { fn pdf(&self, x: &f64) -> f64 { use special_fun::FloatSpecial; let nu = get_nu!(self); let np1o2 = (nu + 1.0) / 2.0; let norm = np1o2.gamma() / (nu * PI).sqrt() / (nu / 2.0).gamma(); norm * (1.0 + x * x / nu).powf(-np1o2) } } impl UnivariateMoments for StudentT { fn mean(&self) -> f64 { if get_nu!(self) <= 1.0 { undefined!("Mean is undefined for nu <= 1."); } 0.0 } fn variance(&self) -> f64 { let nu = get_nu!(self); if nu <= 1.0 { undefined!("Variance is undefined for nu <= 1."); } else if nu <= 2.0 { undefined!("Variance is infinite for 1 < nu <= 2."); } nu / (nu - 2.0) } fn skewness(&self) -> f64 { if get_nu!(self) <= 3.0 { undefined!("Skewness is undefined for nu <= 1."); } 0.0 } fn excess_kurtosis(&self) -> f64 { let nu = get_nu!(self); if nu <= 2.0 { undefined!("Kurtosis is undefined for nu <= 2."); } else if nu <= 4.0 { undefined!("Kurtosis is infinite for 2 < nu <= 4."); } 6.0 / (nu - 4.0) } } impl Quantiles for StudentT { fn quantile(&self, _: Probability) -> f64 { unimplemented!() } fn median(&self) -> f64 { 0.0 } } impl Modes for StudentT { fn modes(&self) -> Vec<f64> { vec![0.0] } } impl Entropy for StudentT { fn entropy(&self) -> f64 { use special_fun::FloatSpecial; let nu = get_nu!(self); let no2 = nu / 2.0; let np1o2 = (nu + 1.0) / 2.0; np1o2 * (np1o2.gamma() - no2.gamma()) + (nu.sqrt() * (no2.beta(0.5))) } } impl fmt::Display for StudentT { fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { write!(f, "StudentT({})", get_nu!(self)) } }