ngboost_rs/dist/
cauchy.rs

1use crate::dist::{Distribution, DistributionMethods, RegressionDistn};
2use crate::scores::{LogScore, Scorable};
3use ndarray::{array, Array1, Array2, Array3};
4use rand::prelude::*;
5use statrs::distribution::{Cauchy as CauchyDist, Continuous, ContinuousCDF};
6
7/// The Cauchy distribution.
8///
9/// The Cauchy distribution is equivalent to the Student's T distribution with df=1.
10/// It has two parameters: loc (median) and log(scale).
11///
12/// Note: The Cauchy distribution has no defined mean or variance.
13/// The `predict` method returns the median (loc).
14#[derive(Debug, Clone)]
15pub struct Cauchy {
16    /// The location parameter (median).
17    pub loc: Array1<f64>,
18    /// The scale parameter.
19    pub scale: Array1<f64>,
20    /// The variance (scale^2) - used in gradient computations.
21    pub var: Array1<f64>,
22    /// The parameters of the distribution, stored as a 2D array.
23    _params: Array2<f64>,
24}
25
26/// Fixed df=1 for Cauchy (T distribution with df=1)
27const CAUCHY_DF: f64 = 1.0;
28
29impl Distribution for Cauchy {
30    fn from_params(params: &Array2<f64>) -> Self {
31        let loc = params.column(0).to_owned();
32        let scale = params.column(1).mapv(f64::exp);
33        let var = &scale * &scale;
34        Cauchy {
35            loc,
36            scale,
37            var,
38            _params: params.clone(),
39        }
40    }
41
42    fn fit(y: &Array1<f64>) -> Array1<f64> {
43        // For Cauchy, use median and interquartile range for robust estimation
44        let n = y.len();
45        if n == 0 {
46            return array![0.0, 0.0];
47        }
48
49        let mut sorted: Vec<f64> = y.to_vec();
50        sorted.sort_by(|a, b| a.partial_cmp(b).unwrap());
51
52        // Median
53        let median = if n % 2 == 0 {
54            (sorted[n / 2 - 1] + sorted[n / 2]) / 2.0
55        } else {
56            sorted[n / 2]
57        };
58
59        // IQR-based scale estimate (half the IQR for Cauchy)
60        let q1_idx = n / 4;
61        let q3_idx = 3 * n / 4;
62        let iqr = sorted[q3_idx] - sorted[q1_idx];
63        let scale = (iqr / 2.0).max(1e-6);
64
65        array![median, scale.ln()]
66    }
67
68    fn n_params(&self) -> usize {
69        2
70    }
71
72    fn predict(&self) -> Array1<f64> {
73        // Return median (Cauchy has no mean)
74        self.loc.clone()
75    }
76
77    fn params(&self) -> &Array2<f64> {
78        &self._params
79    }
80}
81
82impl RegressionDistn for Cauchy {}
83
84impl DistributionMethods for Cauchy {
85    fn mean(&self) -> Array1<f64> {
86        // Cauchy has no defined mean, return NaN
87        Array1::from_elem(self.loc.len(), f64::NAN)
88    }
89
90    fn variance(&self) -> Array1<f64> {
91        // Cauchy has no defined variance, return NaN
92        Array1::from_elem(self.loc.len(), f64::NAN)
93    }
94
95    fn std(&self) -> Array1<f64> {
96        // Cauchy has no defined std, return NaN
97        Array1::from_elem(self.loc.len(), f64::NAN)
98    }
99
100    fn pdf(&self, y: &Array1<f64>) -> Array1<f64> {
101        let mut result = Array1::zeros(y.len());
102        for i in 0..y.len() {
103            if let Ok(d) = CauchyDist::new(self.loc[i], self.scale[i]) {
104                result[i] = d.pdf(y[i]);
105            }
106        }
107        result
108    }
109
110    fn logpdf(&self, y: &Array1<f64>) -> Array1<f64> {
111        let mut result = Array1::zeros(y.len());
112        for i in 0..y.len() {
113            if let Ok(d) = CauchyDist::new(self.loc[i], self.scale[i]) {
114                result[i] = d.ln_pdf(y[i]);
115            }
116        }
117        result
118    }
119
120    fn cdf(&self, y: &Array1<f64>) -> Array1<f64> {
121        let mut result = Array1::zeros(y.len());
122        for i in 0..y.len() {
123            if let Ok(d) = CauchyDist::new(self.loc[i], self.scale[i]) {
124                result[i] = d.cdf(y[i]);
125            }
126        }
127        result
128    }
129
130    fn ppf(&self, q: &Array1<f64>) -> Array1<f64> {
131        // Inverse CDF for Cauchy: loc + scale * tan(pi * (q - 0.5))
132        let mut result = Array1::zeros(q.len());
133        for i in 0..q.len() {
134            let q_clamped = q[i].clamp(1e-15, 1.0 - 1e-15);
135            result[i] =
136                self.loc[i] + self.scale[i] * (std::f64::consts::PI * (q_clamped - 0.5)).tan();
137        }
138        result
139    }
140
141    fn sample(&self, n_samples: usize) -> Array2<f64> {
142        let n_obs = self.loc.len();
143        let mut samples = Array2::zeros((n_samples, n_obs));
144        let mut rng = rand::rng();
145
146        for i in 0..n_obs {
147            for s in 0..n_samples {
148                // Use inverse CDF method: loc + scale * tan(pi * (u - 0.5))
149                let u: f64 = rng.random();
150                samples[[s, i]] =
151                    self.loc[i] + self.scale[i] * (std::f64::consts::PI * (u - 0.5)).tan();
152            }
153        }
154        samples
155    }
156
157    fn median(&self) -> Array1<f64> {
158        // Median of Cauchy is loc
159        self.loc.clone()
160    }
161
162    fn mode(&self) -> Array1<f64> {
163        // Mode of Cauchy is loc
164        self.loc.clone()
165    }
166}
167
168impl Scorable<LogScore> for Cauchy {
169    fn score(&self, y: &Array1<f64>) -> Array1<f64> {
170        let mut scores = Array1::zeros(y.len());
171        for (i, &y_i) in y.iter().enumerate() {
172            let d = CauchyDist::new(self.loc[i], self.scale[i]).unwrap();
173            scores[i] = -d.ln_pdf(y_i);
174        }
175        scores
176    }
177
178    fn d_score(&self, y: &Array1<f64>) -> Array2<f64> {
179        // Same as TFixedDf with df=1
180        let n_obs = y.len();
181        let mut d_params = Array2::zeros((n_obs, 2));
182
183        for i in 0..n_obs {
184            let loc_i = self.loc[i];
185            let var_i = self.var[i];
186            let y_i = y[i];
187
188            let diff = y_i - loc_i;
189            let diff_sq = diff * diff;
190            let denom = CAUCHY_DF * var_i + diff_sq;
191
192            // d/d(loc): -(df + 1) * (y - loc) / (df * var + (y - loc)^2)
193            d_params[[i, 0]] = -(CAUCHY_DF + 1.0) * diff / denom;
194
195            // d/d(log(scale)): 1 - (df + 1) * (y - loc)^2 / (df * var + (y - loc)^2)
196            d_params[[i, 1]] = 1.0 - (CAUCHY_DF + 1.0) * diff_sq / denom;
197        }
198
199        d_params
200    }
201
202    fn metric(&self) -> Array3<f64> {
203        // Fisher Information Matrix for Cauchy (T with df=1)
204        // FI[0, 0] = (df + 1) / ((df + 3) * var) = 2 / (4 * var) = 0.5 / var
205        // FI[1, 1] = df / (2 * (df + 3) * var) = 1 / (8 * var)
206        let n_obs = self.loc.len();
207        let mut fi = Array3::zeros((n_obs, 2, 2));
208
209        for i in 0..n_obs {
210            let var_i = self.var[i];
211            fi[[i, 0, 0]] = (CAUCHY_DF + 1.0) / ((CAUCHY_DF + 3.0) * var_i);
212            fi[[i, 1, 1]] = CAUCHY_DF / (2.0 * (CAUCHY_DF + 3.0) * var_i);
213        }
214
215        fi
216    }
217}
218
219/// The Cauchy distribution with fixed variance=1.
220///
221/// Has one parameter: loc (median).
222#[derive(Debug, Clone)]
223pub struct CauchyFixedVar {
224    /// The location parameter (median).
225    pub loc: Array1<f64>,
226    /// The scale parameter (fixed at 1.0).
227    pub scale: Array1<f64>,
228    /// The variance (fixed at 1.0).
229    pub var: Array1<f64>,
230    /// The parameters of the distribution, stored as a 2D array.
231    _params: Array2<f64>,
232}
233
234impl Distribution for CauchyFixedVar {
235    fn from_params(params: &Array2<f64>) -> Self {
236        let loc = params.column(0).to_owned();
237        let n = loc.len();
238        let scale = Array1::ones(n);
239        let var = Array1::ones(n);
240        CauchyFixedVar {
241            loc,
242            scale,
243            var,
244            _params: params.clone(),
245        }
246    }
247
248    fn fit(y: &Array1<f64>) -> Array1<f64> {
249        // Median estimation
250        let n = y.len();
251        if n == 0 {
252            return array![0.0];
253        }
254
255        let mut sorted: Vec<f64> = y.to_vec();
256        sorted.sort_by(|a, b| a.partial_cmp(b).unwrap());
257
258        let median = if n % 2 == 0 {
259            (sorted[n / 2 - 1] + sorted[n / 2]) / 2.0
260        } else {
261            sorted[n / 2]
262        };
263
264        array![median]
265    }
266
267    fn n_params(&self) -> usize {
268        1
269    }
270
271    fn predict(&self) -> Array1<f64> {
272        self.loc.clone()
273    }
274
275    fn params(&self) -> &Array2<f64> {
276        &self._params
277    }
278}
279
280impl RegressionDistn for CauchyFixedVar {}
281
282impl DistributionMethods for CauchyFixedVar {
283    fn mean(&self) -> Array1<f64> {
284        Array1::from_elem(self.loc.len(), f64::NAN)
285    }
286
287    fn variance(&self) -> Array1<f64> {
288        Array1::from_elem(self.loc.len(), f64::NAN)
289    }
290
291    fn pdf(&self, y: &Array1<f64>) -> Array1<f64> {
292        let mut result = Array1::zeros(y.len());
293        for i in 0..y.len() {
294            if let Ok(d) = CauchyDist::new(self.loc[i], self.scale[i]) {
295                result[i] = d.pdf(y[i]);
296            }
297        }
298        result
299    }
300
301    fn cdf(&self, y: &Array1<f64>) -> Array1<f64> {
302        let mut result = Array1::zeros(y.len());
303        for i in 0..y.len() {
304            if let Ok(d) = CauchyDist::new(self.loc[i], self.scale[i]) {
305                result[i] = d.cdf(y[i]);
306            }
307        }
308        result
309    }
310
311    fn ppf(&self, q: &Array1<f64>) -> Array1<f64> {
312        let mut result = Array1::zeros(q.len());
313        for i in 0..q.len() {
314            let q_clamped = q[i].clamp(1e-15, 1.0 - 1e-15);
315            result[i] =
316                self.loc[i] + self.scale[i] * (std::f64::consts::PI * (q_clamped - 0.5)).tan();
317        }
318        result
319    }
320
321    fn sample(&self, n_samples: usize) -> Array2<f64> {
322        let n_obs = self.loc.len();
323        let mut samples = Array2::zeros((n_samples, n_obs));
324        let mut rng = rand::rng();
325
326        for i in 0..n_obs {
327            for s in 0..n_samples {
328                let u: f64 = rng.random();
329                samples[[s, i]] =
330                    self.loc[i] + self.scale[i] * (std::f64::consts::PI * (u - 0.5)).tan();
331            }
332        }
333        samples
334    }
335
336    fn median(&self) -> Array1<f64> {
337        self.loc.clone()
338    }
339
340    fn mode(&self) -> Array1<f64> {
341        self.loc.clone()
342    }
343}
344
345impl Scorable<LogScore> for CauchyFixedVar {
346    fn score(&self, y: &Array1<f64>) -> Array1<f64> {
347        let mut scores = Array1::zeros(y.len());
348        for (i, &y_i) in y.iter().enumerate() {
349            let d = CauchyDist::new(self.loc[i], self.scale[i]).unwrap();
350            scores[i] = -d.ln_pdf(y_i);
351        }
352        scores
353    }
354
355    fn d_score(&self, y: &Array1<f64>) -> Array2<f64> {
356        let n_obs = y.len();
357        let mut d_params = Array2::zeros((n_obs, 1));
358
359        for i in 0..n_obs {
360            let loc_i = self.loc[i];
361            let var_i = self.var[i];
362            let y_i = y[i];
363
364            let diff = y_i - loc_i;
365            let diff_sq = diff * diff;
366
367            // d/d(loc) for fixed var case
368            let num = (CAUCHY_DF + 1.0) * (2.0 / (CAUCHY_DF * var_i)) * diff;
369            let den = 2.0 * (1.0 + (1.0 / (CAUCHY_DF * var_i)) * diff_sq);
370            d_params[[i, 0]] = -num / den;
371        }
372
373        d_params
374    }
375
376    fn metric(&self) -> Array3<f64> {
377        let n_obs = self.loc.len();
378        let mut fi = Array3::zeros((n_obs, 1, 1));
379
380        for i in 0..n_obs {
381            let var_i = self.var[i];
382            fi[[i, 0, 0]] = (CAUCHY_DF + 1.0) / ((CAUCHY_DF + 3.0) * var_i);
383        }
384
385        fi
386    }
387}
388
389#[cfg(test)]
390mod tests {
391    use super::*;
392    use approx::assert_relative_eq;
393
394    #[test]
395    fn test_cauchy_distribution_methods() {
396        let params = Array2::from_shape_vec((2, 2), vec![0.0, 0.0, 5.0, 1.0_f64.ln()]).unwrap();
397        let dist = Cauchy::from_params(&params);
398
399        // Mean and variance should be NaN for Cauchy
400        assert!(dist.mean()[0].is_nan());
401        assert!(dist.variance()[0].is_nan());
402
403        // Median should be loc
404        let median = dist.median();
405        assert_relative_eq!(median[0], 0.0, epsilon = 1e-10);
406        assert_relative_eq!(median[1], 5.0, epsilon = 1e-10);
407
408        // Mode should be loc
409        let mode = dist.mode();
410        assert_relative_eq!(mode[0], 0.0, epsilon = 1e-10);
411        assert_relative_eq!(mode[1], 5.0, epsilon = 1e-10);
412    }
413
414    #[test]
415    fn test_cauchy_cdf_ppf() {
416        // Create 3 observations for ppf inverse test
417        let params = Array2::from_shape_vec((3, 2), vec![0.0, 0.0, 0.0, 0.0, 0.0, 0.0]).unwrap();
418        let dist = Cauchy::from_params(&params);
419
420        // CDF at loc should be 0.5
421        let y = Array1::from_vec(vec![0.0, 0.0, 0.0]);
422        let cdf = dist.cdf(&y);
423        assert_relative_eq!(cdf[0], 0.5, epsilon = 1e-10);
424
425        // PPF at 0.5 should be loc
426        let q = Array1::from_vec(vec![0.5, 0.5, 0.5]);
427        let ppf = dist.ppf(&q);
428        assert_relative_eq!(ppf[0], 0.0, epsilon = 1e-10);
429
430        // PPF inverse test
431        let q = Array1::from_vec(vec![0.25, 0.5, 0.75]);
432        let ppf = dist.ppf(&q);
433        let cdf_of_ppf = dist.cdf(&ppf);
434        assert_relative_eq!(cdf_of_ppf[0], 0.25, epsilon = 1e-6);
435        assert_relative_eq!(cdf_of_ppf[1], 0.5, epsilon = 1e-6);
436        assert_relative_eq!(cdf_of_ppf[2], 0.75, epsilon = 1e-6);
437    }
438
439    #[test]
440    fn test_cauchy_pdf() {
441        let params = Array2::from_shape_vec((1, 2), vec![0.0, 0.0]).unwrap();
442        let dist = Cauchy::from_params(&params);
443
444        // PDF at loc should be 1/(pi*scale) = 1/pi for scale=1
445        let y = Array1::from_vec(vec![0.0]);
446        let pdf = dist.pdf(&y);
447        assert_relative_eq!(pdf[0], 1.0 / std::f64::consts::PI, epsilon = 1e-10);
448    }
449
450    #[test]
451    fn test_cauchy_sample() {
452        let params = Array2::from_shape_vec((1, 2), vec![5.0, 0.0]).unwrap();
453        let dist = Cauchy::from_params(&params);
454
455        let samples = dist.sample(1000);
456        assert_eq!(samples.shape(), &[1000, 1]);
457
458        // Check that sample median is close to loc = 5
459        // (mean is not defined for Cauchy, so we check median)
460        let mut sample_vec: Vec<f64> = samples.column(0).to_vec();
461        sample_vec.sort_by(|a, b| a.partial_cmp(b).unwrap());
462        let sample_median = sample_vec[500];
463        assert!((sample_median - 5.0).abs() < 1.0);
464    }
465
466    #[test]
467    fn test_cauchy_fit() {
468        let y = Array1::from_vec(vec![1.0, 2.0, 3.0, 4.0, 5.0]);
469        let params = Cauchy::fit(&y);
470        assert_eq!(params.len(), 2);
471        // Median should be 3.0
472        assert_relative_eq!(params[0], 3.0, epsilon = 1e-10);
473    }
474
475    #[test]
476    fn test_cauchy_logscore() {
477        let params = Array2::from_shape_vec((1, 2), vec![0.0, 0.0]).unwrap();
478        let dist = Cauchy::from_params(&params);
479
480        let y = Array1::from_vec(vec![0.0]);
481        let score = Scorable::<LogScore>::score(&dist, &y);
482
483        // Score at loc should be -ln(1/(pi*scale)) = ln(pi) for scale=1
484        assert_relative_eq!(score[0], std::f64::consts::PI.ln(), epsilon = 1e-10);
485    }
486
487    #[test]
488    fn test_cauchy_fixed_var() {
489        let params = Array2::from_shape_vec((1, 1), vec![5.0]).unwrap();
490        let dist = CauchyFixedVar::from_params(&params);
491
492        assert_relative_eq!(dist.median()[0], 5.0, epsilon = 1e-10);
493        assert_relative_eq!(dist.scale[0], 1.0, epsilon = 1e-10);
494
495        let y = Array1::from_vec(vec![5.0]);
496        let cdf = dist.cdf(&y);
497        assert_relative_eq!(cdf[0], 0.5, epsilon = 1e-10);
498    }
499
500    #[test]
501    fn test_cauchy_interval() {
502        let params = Array2::from_shape_vec((1, 2), vec![0.0, 0.0]).unwrap();
503        let dist = Cauchy::from_params(&params);
504
505        let (lower, upper) = dist.interval(0.5);
506        // For Cauchy with loc=0, scale=1, the 25% and 75% quantiles are -1 and 1
507        assert_relative_eq!(lower[0], -1.0, epsilon = 1e-10);
508        assert_relative_eq!(upper[0], 1.0, epsilon = 1e-10);
509    }
510}
ngboost_rs/dist/cauchy.rs

ngboost_rs/dist/
cauchy.rs
