# Statistics

Ten probability distributions, each implementing `ContinuousDistribution<T>` or `DiscreteDistribution<T>` traits.

Requires the `stats` Cargo feature (implies `special`):

```toml
numeris = { version = "0.5", features = ["stats"] }
```

All distributions work with `f32` and `f64`, are no-std compatible, and have no heap allocation.

## Distribution Traits

```rust
pub trait ContinuousDistribution<T> {
    fn pdf(&self, x: T) -> T;       // probability density function
    fn ln_pdf(&self, x: T) -> T;    // log-PDF (avoids underflow)
    fn cdf(&self, x: T) -> T;       // cumulative distribution function
    fn quantile(&self, p: T) -> T;  // inverse CDF: x such that P(X ≤ x) = p
    fn mean(&self) -> T;
    fn variance(&self) -> T;
}

pub trait DiscreteDistribution<T> {
    fn pmf(&self, k: u64) -> T;     // probability mass function P(X = k)
    fn ln_pmf(&self, k: u64) -> T;  // log-PMF (avoids underflow)
    fn cdf(&self, k: u64) -> T;     // cumulative P(X ≤ k)
    fn quantile(&self, p: T) -> u64; // smallest k such that P(X ≤ k) ≥ p
    fn mean(&self) -> T;
    fn variance(&self) -> T;
}
```

## Continuous Distributions

### Normal (Gaussian)

`Normal<T> { mean: T, std_dev: T }`

$$
f(x) = \frac{1}{\sigma\sqrt{2\pi}}\,
        \exp\!\left(-\frac{(x-\mu)^2}{2\sigma^2}\right),
\qquad x \in \mathbb{R}.
$$

Mean $\mu$, variance $\sigma^2$.

![Normal distribution PDF plot](includes/plot_normal_pdf.svg)

```rust
use numeris::stats::Normal;

let n = Normal::new(0.0_f64, 1.0).unwrap();  // standard normal N(0,1)
let p = n.pdf(0.0);      // 1/√(2π) ≈ 0.3989
let c = n.cdf(1.645);    // ≈ 0.95 (95th percentile)

let n2 = Normal::new(100.0_f64, 15.0).unwrap(); // IQ distribution
let prob_above_130 = 1.0 - n2.cdf(130.0);       // P(IQ > 130) ≈ 2.3%
```

### Uniform

`Uniform<T> { a: T, b: T }` — continuous uniform on [a, b].

$$
f(x) = \frac{1}{b - a} \quad \text{for } x \in [a, b], \qquad 0 \text{ otherwise.}
$$

Mean $\dfrac{a+b}{2}$, variance $\dfrac{(b-a)^2}{12}$.

```rust
use numeris::stats::Uniform;

let u = Uniform::new(0.0_f64, 1.0).unwrap();
let p = u.pdf(0.5);      // 1.0 (constant density)
let c = u.cdf(0.3);      // 0.3
assert_eq!(u.mean(), 0.5);
assert!((u.variance() - 1.0/12.0).abs() < 1e-10);
```

### Exponential

`Exponential<T> { rate: T }` — rate λ (mean = 1/λ).

$$
f(x) = \lambda e^{-\lambda x}, \qquad x \ge 0.
$$

Mean $1/\lambda$, variance $1/\lambda^2$.

```rust
use numeris::stats::Exponential;

let e = Exponential::new(2.0_f64).unwrap();  // λ=2, mean=0.5
let p = e.pdf(1.0);      // 2·e⁻² ≈ 0.2707
let c = e.cdf(1.0);      // 1 - e⁻² ≈ 0.8647
assert!((e.mean() - 0.5).abs() < 1e-12);
```

### Gamma

`Gamma<T> { shape: T, rate: T }` — shape α, rate β (mean = α/β).

$$
f(x) = \frac{\beta^{\alpha}}{\Gamma(\alpha)}\,
        x^{\alpha-1} e^{-\beta x}, \qquad x > 0.
$$

Mean $\alpha/\beta$, variance $\alpha/\beta^2$. The chi-squared and exponential
distributions are special cases.

![Gamma distribution PDF plot](includes/plot_gamma_pdf.svg)

```rust
use numeris::stats::Gamma;

let g = Gamma::new(2.0_f64, 1.0).unwrap();  // Gamma(α=2, β=1)
let p = g.pdf(1.0);   // 1·e⁻¹ ≈ 0.3679
let c = g.cdf(2.0);   // uses regularized incomplete gamma P(2, 2)
assert!((g.mean() - 2.0).abs() < 1e-10);
assert!((g.variance() - 2.0).abs() < 1e-10);
```

### Beta

`Beta<T> { alpha: T, beta: T }` — shape parameters α, β. Support: [0, 1].

$$
f(x) = \frac{x^{\alpha-1}(1-x)^{\beta-1}}{B(\alpha, \beta)},
\qquad x \in [0, 1], \quad
B(\alpha, \beta) = \frac{\Gamma(\alpha)\,\Gamma(\beta)}{\Gamma(\alpha+\beta)}.
$$

Mean $\dfrac{\alpha}{\alpha+\beta}$, variance
$\dfrac{\alpha\beta}{(\alpha+\beta)^2(\alpha+\beta+1)}$.

![Beta distribution PDF plot](includes/plot_beta_pdf.svg)

```rust
use numeris::stats::Beta;

let b = Beta::new(2.0_f64, 5.0).unwrap();
let p = b.pdf(0.3);    // uses regularized incomplete beta
let c = b.cdf(0.5);
assert!((b.mean() - 2.0/7.0).abs() < 1e-10);
```

### Chi-Squared

`ChiSquared<T> { k: T }` — k degrees of freedom. Special case of Gamma(k/2, 1/2).

$$
f(x) = \frac{1}{2^{k/2}\,\Gamma(k/2)}\,
        x^{k/2 - 1} e^{-x/2}, \qquad x > 0.
$$

Mean $k$, variance $2k$.

```rust
use numeris::stats::ChiSquared;

let chi = ChiSquared::new(3.0_f64).unwrap();  // χ²(3)
let p_val = 1.0 - chi.cdf(7.815);    // p-value for χ²=7.815 with 3 df ≈ 0.05
assert_eq!(chi.mean(), 3.0);
assert_eq!(chi.variance(), 6.0);
```

### Student's t

`StudentT<T> { nu: T }` — ν degrees of freedom.

$$
f(x) = \frac{\Gamma\!\left(\frac{\nu+1}{2}\right)}
            {\sqrt{\nu\pi}\;\Gamma\!\left(\frac{\nu}{2}\right)}
        \left(1 + \frac{x^2}{\nu}\right)^{-\frac{\nu+1}{2}}.
$$

Mean $0$ (for $\nu > 1$), variance $\dfrac{\nu}{\nu - 2}$ (for $\nu > 2$).

```rust
use numeris::stats::StudentT;

let t = StudentT::new(10.0_f64).unwrap();  // t(10)
let c = t.cdf(2.228);   // ≈ 0.975 (95% two-sided CI critical value)
let p = t.pdf(0.0);     // peak of distribution
assert_eq!(t.mean(), 0.0);
assert!((t.variance() - 10.0/8.0).abs() < 1e-10);  // ν/(ν-2)
```

### CDF Comparison

![Continuous CDF comparison plot](includes/plot_continuous_cdf.svg)

## Discrete Distributions

### Bernoulli

`Bernoulli<T> { p: T }` — single trial with success probability p.

$$
P(X = k) = p^k (1-p)^{1-k}, \qquad k \in \{0, 1\}.
$$

Mean $p$, variance $p(1-p)$.

```rust
use numeris::stats::Bernoulli;

let b = Bernoulli::new(0.3_f64).unwrap();
assert!((b.pmf(1) - 0.3).abs() < 1e-12);
assert!((b.pmf(0) - 0.7).abs() < 1e-12);
assert!((b.mean() - 0.3).abs() < 1e-12);
assert!((b.variance() - 0.21).abs() < 1e-12);
```

### Binomial

`Binomial<T> { n: u64, p: T }` — n independent Bernoulli trials.

$$
P(X = k) = \binom{n}{k} p^k (1-p)^{n-k},
\qquad k = 0, 1, \dots, n.
$$

Mean $np$, variance $np(1-p)$.

![Binomial distribution PMF plot](includes/plot_binomial_pmf.svg)

```rust
use numeris::stats::Binomial;

let b = Binomial::new(10, 0.5_f64).unwrap();    // B(10, 0.5)
let p5 = b.pmf(5);    // P(X=5) ≈ 0.2461
let c7  = b.cdf(7);   // P(X≤7) ≈ 0.9453
assert!((b.mean() - 5.0).abs() < 1e-10);
assert!((b.variance() - 2.5).abs() < 1e-10);
```

### Poisson

`Poisson<T> { lambda: T }` — number of events in a fixed interval.

$$
P(X = k) = \frac{\lambda^k e^{-\lambda}}{k!},
\qquad k = 0, 1, 2, \dots
$$

Mean $\lambda$, variance $\lambda$.

![Poisson distribution PMF plot](includes/plot_poisson_pmf.svg)

```rust
use numeris::stats::Poisson;

let pois = Poisson::new(3.0_f64).unwrap();  // Poisson(λ=3)
let p3 = pois.pmf(3);    // P(X=3) = e⁻³·3³/3! ≈ 0.2240
let c5 = pois.cdf(5);    // P(X≤5) ≈ 0.9161
assert_eq!(pois.mean(), 3.0);
assert_eq!(pois.variance(), 3.0);
```

## Distribution Summary Table

| Distribution | Type | Parameters | Mean | Variance |
|---|---|---|---|---|
| `Normal` | Continuous | μ, σ | μ | σ² |
| `Uniform` | Continuous | a, b | (a+b)/2 | (b-a)²/12 |
| `Exponential` | Continuous | λ | 1/λ | 1/λ² |
| `Gamma` | Continuous | α, β | α/β | α/β² |
| `Beta` | Continuous | α, β | α/(α+β) | αβ/((α+β)²(α+β+1)) |
| `ChiSquared` | Continuous | k | k | 2k |
| `StudentT` | Continuous | ν | 0 (ν>1) | ν/(ν-2) (ν>2) |
| `Bernoulli` | Discrete | p | p | p(1-p) |
| `Binomial` | Discrete | n, p | np | np(1-p) |
| `Poisson` | Discrete | λ | λ | λ |

## Random Sampling

Every distribution provides `sample()` and `sample_array()` methods using the built-in [`Rng`] (xoshiro256++ PRNG). No external crate dependency.

```rust
use numeris::stats::{Rng, Normal, Uniform, Poisson};
use numeris::stats::{ContinuousDistribution, DiscreteDistribution};

let mut rng = Rng::new(42);  // seed with any u64

// Sample from a normal distribution
let n = Normal::new(0.0_f64, 1.0).unwrap();
let x = n.sample(&mut rng);

// Fill an array with 100 samples
let samples: [f64; 100] = n.sample_array(&mut rng);

// Works with all distributions
let u = Uniform::new(0.0_f64, 10.0).unwrap();
let y = u.sample(&mut rng);

let pois = Poisson::new(3.0_f64).unwrap();
let k = pois.sample(&mut rng);  // returns T (float), cast to integer if needed
```

The `Rng` struct is deterministic — same seed produces the same sequence, useful for reproducible simulations.

```rust
use numeris::stats::Rng;

let mut rng1 = Rng::new(123);
let mut rng2 = Rng::new(123);
assert_eq!(rng1.next_u64(), rng2.next_u64());
assert_eq!(rng1.next_f64(), rng2.next_f64());
```

### Rng Methods

| Method | Returns |
|---|---|
| `next_u64()` | Uniform `u64` |
| `next_f64()` | Uniform `f64` in [0, 1) |
| `next_f32()` | Uniform `f32` in [0, 1) |

### Sampling Algorithms

| Distribution | Algorithm |
|---|---|
| Normal | Box-Muller transform |
| Uniform | Inverse CDF (affine transform) |
| Exponential | Inverse CDF (−ln(U)/λ) |
| Gamma | Marsaglia & Tsang's method |
| Beta | Ratio of two Gamma variates |
| Chi-Squared | Gamma(k/2, 1/2) |
| Student's t | Normal / √(χ²/ν) |
| Bernoulli | Threshold on uniform |
| Binomial | Sum of Bernoulli (small n) or normal approximation (large n) |
| Poisson | Knuth's method (small λ) or normal approximation (large λ) |

## Quantile (Inverse CDF)

All distributions implement `quantile()` for computing percentiles:

```rust
use numeris::stats::{Normal, ContinuousDistribution};

let n = Normal::new(0.0_f64, 1.0).unwrap();
let z_95 = n.quantile(0.95);  // ≈ 1.645
let z_99 = n.quantile(0.99);  // ≈ 2.326
```

```rust
use numeris::stats::{Binomial, DiscreteDistribution};

let b = Binomial::new(20, 0.5_f64).unwrap();
let k = b.quantile(0.95);  // smallest k with P(X ≤ k) ≥ 0.95
assert_eq!(k, 14);
```

## Error Handling

```rust
use numeris::stats::StatsError;

match Normal::new(0.0_f64, -1.0) {  // std_dev must be > 0
    Err(StatsError::InvalidParameter) => { /* bad parameter */ }
    Ok(n) => { /* use n */ }
}
```
