Add Pareto (Type I) distribution

CHANGELOG.md

@@ -5,6 +5,15 @@ All notable changes to this crate are documented here. The format is based on
 [Semantic Versioning](https://semver.org/spec/v2.0.0.html). Pre-`1.0`, minor versions may
 contain breaking changes.
 
+## [Unreleased]
+
+### Added
+- **Pareto (Type I) distribution.** `ParetoDistribution` (in `distributions::positive`) implements
+  `Pdf`, `Cdf`, `Quantile`, `Moments`, and `Sample`, with scipy-validated golden fixtures
+  (`dist_pareto.json`) and full equivalence, moment-convergence, tail-stress, and parameter-sweep
+  coverage. The mean is defined for `α > 1` and the variance for `α > 2`, reported as `None`
+  (matching scipy's NaN) otherwise.
+
 ## [0.2.1] — 2026-08-02
 
 A release-engineering fix. The library's public API, numerics, and behaviour are identical

README.md

@@ -75,8 +75,8 @@ assert_eq!(z.variance(), Some(1.0));
 ## What's in the box
 
 **Distributions** — `pdf`/`pmf`, `cdf`, `quantile`, `moments`, log-space `cdf`/`sf`, and
-seeded sampling for 14 families: normal, uniform, Cauchy, Laplace, exponential, gamma, beta,
-Weibull, log-normal, chi-squared, Student's t, F, binomial, and Poisson.
+seeded sampling for 15 families: normal, uniform, Cauchy, Laplace, exponential, gamma, beta,
+Weibull, log-normal, Pareto, chi-squared, Student's t, F, binomial, and Poisson.
 
 **Hypothesis tests** — parametric (one-sample / independent / paired / Welch t-tests, ANOVA,
 variance), non-parametric (Mann–Whitney, Wilcoxon, Kruskal–Wallis, Friedman), categorical

reference/gen/gen_distributions.py

@@ -196,6 +196,20 @@ def gen_beta():
     )
 
 
+def gen_pareto():
+    alpha, x_m = 3.0, 1.5
+    d = stats.pareto(b=alpha, scale=x_m)
+    # Support is [x_m, inf); span the body and a stretch of the right tail.
+    xs = np.linspace(x_m, 10.0 * x_m, 41)
+    _continuous(
+        "dist_pareto",
+        d,
+        {"shape_parameter": alpha, "scale_parameter": x_m},
+        xs,
+        library="scipy.stats.pareto",
+    )
+
+
 def gen_chi_squared():
     df = 5
     d = stats.chi2(df)
@@ -272,6 +286,7 @@ def main():
     gen_lognormal()
     gen_gamma()
     gen_beta()
+    gen_pareto()
     gen_chi_squared()
     gen_students_t()
     gen_f()

reference/golden/dist_pareto.json

@@ -0,0 +1,164 @@
+{
+  "_provenance": {
+    "library": "scipy.stats.pareto",
+    "seed": null,
+    "version": "1.17.1"
+  },
+  "cdf": [
+    0.0,
+    0.45600897585189837,
+    0.6719832711468284,
+    0.7872078679890812,
+    0.8542061525003645,
+    0.8957866883777732,
+    0.9229457827263708,
+    0.9414309338013978,
+    0.9544460641399417,
+    0.9638736684765582,
+    0.9708693673190715,
+    0.9761693672855308,
+    0.9802578327048743,
+    0.9834620750496202,
+    0.9860087759952569,
+    0.9880583090379009,
+    0.9897263088682502,
+    0.991097580670177,
+    0.9922352788165788,
+    0.9931870906865294,
+    0.9939894815927873,
+    0.9946706587873467,
+    0.9952526734865969,
+    0.9957529329813488,
+    0.996185302734375,
+    0.9965609194166997,
+    0.996888797454103,
+    0.9971762862578633,
+    0.9974294182516446,
+    0.9976531761821877,
+    0.997851700177906,
+    0.9980284494133543,
+    0.9981863292755474,
+    0.998327792097945,
+    0.998454917485834,
+    0.9985694767693645,
+    0.9986729850272901,
+    0.998766743314212,
+    0.9988518731187762,
+    0.9989293446245765,
+    0.999
+  ],
+  "mean": 2.25,
+  "p": [
+    1e-06,
+    0.01,
+    0.1,
+    0.25,
+    0.5,
+    0.75,
+    0.9,
+    0.99,
+    0.999999
+  ],
+  "params": {
+    "shape_parameter": 3.0,
+    "scale_parameter": 1.5
+  },
+  "pdf": [
+    2.0,
+    0.8881486108540435,
+    0.452436867383685,
+    0.2540801576249777,
+    0.15346720789435317,
+    0.09808311682091929,
+    0.06557805725415249,
+    0.045490536853283275,
+    0.03253852561432736,
+    0.02388517786673836,
+    0.017926543188263717,
+    0.013715472065881537,
+    0.010671441781149007,
+    0.008426968127582046,
+    0.006742758556502664,
+    0.005459058725531029,
+    0.00446682223119557,
+    0.003690122001999209,
+    0.0030751371023450117,
+    0.0025830935785670452,
+    0.002185643057168226,
+    0.0018617785895732007,
+    0.0015957400045052722,
+    0.0013755682651501848,
+    0.0011920928955078125,
+    0.0010382130062793384,
+    0.0009083803053713994,
+    0.0007982229659750625,
+    0.0007042689721521838,
+    0.0006237405495846521,
+    0.0005543999540887538,
+    0.0004944327490020461,
+    0.0004423587132811302,
+    0.0003969633001911029,
+    0.00035724451194589426,
+    0.0003223714322558984,
+    0.00029165164235383227,
+    0.00026450545539688795,
+    0.0002404454201515884,
+    0.00021905992336031734,
+    0.0002
+  ],
+  "ppf": [
+    1.5000005000003334,
+    1.5050335947719786,
+    1.5536162529769295,
+    1.6509636244473134,
+    1.8898815748423097,
+    2.381101577952299,
+    3.231652035047826,
+    6.962383250419166,
+    149.9999999985622
+  ],
+  "variance": 1.6875,
+  "x": [
+    1.5,
+    1.8375,
+    2.175,
+    2.5125,
+    2.85,
+    3.1875,
+    3.525,
+    3.8625,
+    4.2,
+    4.5375,
+    4.875,
+    5.2125,
+    5.55,
+    5.8875,
+    6.225,
+    6.5625,
+    6.9,
+    7.2375,
+    7.575,
+    7.9125,
+    8.25,
+    8.5875,
+    8.925,
+    9.2625,
+    9.6,
+    9.9375,
+    10.275,
+    10.6125,
+    10.95,
+    11.2875,
+    11.625,
+    11.9625,
+    12.3,
+    12.6375,
+    12.975,
+    13.3125,
+    13.65,
+    13.9875,
+    14.325,
+    14.6625,
+    15.0
+  ]
+}

src/distributions/positive/mod.rs

@@ -6,4 +6,5 @@ pub mod beta;
 pub mod exponential;
 pub mod gamma;
 pub mod lognormal;
+pub mod pareto;
 pub mod weibull;

src/distributions/positive/pareto.rs

@@ -0,0 +1,152 @@
+//! Pareto (Type I) distribution numerics, for the [`ParetoDistribution`].
+//!
+//! Equivalent to `scipy.stats.pareto(b = shape_parameter, scale = scale_parameter)`:
+//! a power-law distribution on `[scale_parameter, ∞)` whose density falls off as
+//! `x ** -(α + 1)` above the scale. Every function is closed form; sampling is the
+//! inverse-CDF transform of a single uniform draw.
+//!
+//! The mean is defined only for `α > 1` (`α · x_m / (α - 1)`) and the variance only
+//! for `α > 2` (`α · x_m² / ((α - 1)² · (α - 2))`); both report `None` (matching
+//! scipy's NaN) otherwise.
+//!
+//! # Examples
+//!
+//! ```
+//! use stats_claw::distributions::{Cdf, Moments, Pdf};
+//! use stats_claw::distributions::ParetoDistribution;
+//!
+//! let d = ParetoDistribution { shape_parameter: 2.0, scale_parameter: 1.0, ..Default::default() };
+//! // pdf at the scale is the shape.
+//! assert!((d.pdf(1.0) - 2.0).abs() < 1e-12, "pdf was {}", d.pdf(1.0));
+//! // cdf at 2 is 1 - (1/2)² = 0.75.
+//! assert!((d.cdf(2.0) - 0.75).abs() < 1e-12);
+//! // mean is α·x_m/(α-1) = 2.
+//! assert_eq!(d.mean(), Some(2.0));
+//! ```
+
+use super::super::{Cdf, Moments, Pdf, Quantile, Sample};
+use crate::distributions::ParetoDistribution;
+use crate::rng::SplitMix64;
+
+impl Pdf for ParetoDistribution {
+    fn pdf(&self, x: f64) -> f64 {
+        if x < self.scale_parameter {
+            0.0
+        } else {
+            let alpha = self.shape_parameter;
+            let x_m = self.scale_parameter;
+            alpha * x_m.powf(alpha) / x.powf(alpha + 1.0)
+        }
+    }
+}
+
+impl Cdf for ParetoDistribution {
+    fn cdf(&self, x: f64) -> f64 {
+        if x < self.scale_parameter {
+            0.0
+        } else {
+            let alpha = self.shape_parameter;
+            let x_m = self.scale_parameter;
+            // `1 - (x_m / x) ** alpha` is exact; no `1 - cdf` cancellation because we
+            // evaluate the survival directly, matching scipy's closed form.
+            1.0 - (x_m / x).powf(alpha)
+        }
+    }
+}
+
+impl Quantile for ParetoDistribution {
+    fn quantile(&self, p: f64) -> f64 {
+        let alpha = self.shape_parameter;
+        let x_m = self.scale_parameter;
+        // Direct inverse CDF: x_m / (1 - p) ** (1 / alpha).
+        x_m / (1.0 - p).powf(1.0 / alpha)
+    }
+}
+
+impl Moments for ParetoDistribution {
+    fn mean(&self) -> Option<f64> {
+        let alpha = self.shape_parameter;
+        let x_m = self.scale_parameter;
+        if alpha > 1.0 {
+            Some(alpha * x_m / (alpha - 1.0))
+        } else {
+            None
+        }
+    }
+
+    fn variance(&self) -> Option<f64> {
+        let alpha = self.shape_parameter;
+        let x_m = self.scale_parameter;
+        if alpha > 2.0 {
+            Some(alpha * x_m * x_m / ((alpha - 1.0) * (alpha - 1.0) * (alpha - 2.0)))
+        } else {
+            None
+        }
+    }
+}
+
+impl Sample for ParetoDistribution {
+    fn sample(&self, rng: &mut SplitMix64) -> f64 {
+        self.quantile(rng.next_f64())
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    /// The density at the scale equals the shape parameter.
+    #[test]
+    fn density_at_scale_is_shape() {
+        let d = ParetoDistribution {
+            shape_parameter: 2.0,
+            scale_parameter: 1.0,
+            ..Default::default()
+        };
+        assert!((d.pdf(1.0) - 2.0).abs() < 1e-12, "was {}", d.pdf(1.0));
+    }
+
+    /// The CDF at `2·x_m` with `α = 2` is `1 - (1/2)² = 0.75`.
+    #[test]
+    fn cdf_at_double_scale() {
+        let d = ParetoDistribution {
+            shape_parameter: 2.0,
+            scale_parameter: 1.0,
+            ..Default::default()
+        };
+        assert!((d.cdf(2.0) - 0.75).abs() < 1e-12);
+    }
+
+    /// The mean is `α·x_m/(α-1)` for `α > 1`.
+    #[test]
+    fn mean_is_defined_above_one() {
+        let d = ParetoDistribution {
+            shape_parameter: 3.0,
+            scale_parameter: 1.5,
+            ..Default::default()
+        };
+        assert_eq!(d.mean(), Some(2.25));
+    }
+
+    /// The mean is undefined for `α ≤ 1`, matching scipy's NaN.
+    #[test]
+    fn mean_undefined_at_or_below_one() {
+        let d = ParetoDistribution {
+            shape_parameter: 1.0,
+            scale_parameter: 1.0,
+            ..Default::default()
+        };
+        assert_eq!(d.mean(), None);
+    }
+
+    /// The variance is undefined for `α ≤ 2`.
+    #[test]
+    fn variance_undefined_at_or_below_two() {
+        let d = ParetoDistribution {
+            shape_parameter: 2.0,
+            scale_parameter: 1.0,
+            ..Default::default()
+        };
+        assert_eq!(d.variance(), None);
+    }
+}

src/distributions/types.rs

@@ -190,3 +190,16 @@ pub struct WeibullDistribution {
     /// Free-text description.
     pub description: String,
 }
+
+/// Pareto (Type I) power-law distribution on `[scale_parameter, ∞)`.
+#[derive(Debug, Clone, Default)]
+pub struct ParetoDistribution {
+    /// Shape parameter `α > 0`; controls the tail index.
+    pub shape_parameter: f64,
+    /// Scale (minimum) parameter `x_m > 0`; the lower bound of the support.
+    pub scale_parameter: f64,
+    /// Unique name identifying a distribution.
+    pub distribution_name: String,
+    /// Free-text description.
+    pub description: String,
+}