> For the complete documentation index, see [llms.txt](https://docs.algenta.ai/llms.txt). Markdown versions of documentation pages are available by appending `.md` to page URLs; this page is available as [Markdown](https://docs.algenta.ai/numerical-and-mathematical-computing/probability-distributions-and-dependence-kernels.md).

# Probability distributions & dependence

Sampling and MLE/method-of-moments fitting for twenty-six probability distributions, core probability theory (Bayes, combinatorics, entropy), a Gaussian copula with rank-correlation measures, and Chol

Five modules for the mechanics underneath a Monte Carlo model: what shape a random draw takes (`distributions`), how to recover that shape's parameters from observed data (`distribution_fit`), what a probability calculation is logically allowed to say (`probability`), and how two or more variables move together when a single correlation number does not capture the whole relationship (`copula`, `correlation`). None of the five imports another — each is a flat set of `Float64` formulas over `std.math`/`std.random`, with no shared struct or status-code vocabulary.

**`correlation` is not the module that measures correlation.** Reading its source is what catches this: `correlation.mojo`'s docstring is "Cholesky decomposition for correlated Normal sampling + inverse CDF for QMC/LHS," and every function in it either factors a correlation matrix (`cholesky_lower`) so it can be used to *generate* correlated draws (`sample_correlated_normal`), or maps a uniform variate to a distribution's quantile (`probit`, `icdf`) for stratified/QMC sampling. It has no Pearson, Spearman, or Kendall calculation anywhere in it. Those live one file over, in `copula.mojo`, alongside the Gaussian copula sampler — `rank_correlation` (Spearman) and `kendall_tau` sit in the module named for dependence *modeling*, not the one named for correlation itself.

```python
from algenta import MojoRuntime

runtime = MojoRuntime(mode="local")

# A screening test: 99% sensitive, 5% false-positive rate, 1% base rate in the population.
# P(positive) = 0.99*0.01 + 0.05*0.99 = 0.0594
runtime.execute("probability", "bayes", [0.01, 0.99, 0.0594])["result"]
# 0.16666666666666669 -- a positive result is still only 1-in-6 odds of actually having it

# Which distribution family (by AIC) best explains [1, 2, 3, 4, 5]?
runtime.execute("distribution_fit", "best_fit", [[1.0, 2.0, 3.0, 4.0, 5.0]])["result"]
# 2 -- uniform, not normal (0) or exponential (1); see the worked comparison below
```

## The modules

| Module             | What it does                                                                                                                                                                                                                                                                                    |
| ------------------ | ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- |
| `distributions`    | Twenty-six single-draw Monte Carlo samplers (normal through hypergeometric), a universal `sample_one` dispatcher plus a `sample_one_simd4` batch variant, and SIMD4 batch samplers for the ten distributions used most often in hot loops — all drawing from the process's shared, unseeded RNG |
| `distribution_fit` | MLE fitting for normal and exponential, min/max fitting for uniform, a method-of-moments approximation for Weibull, their log-likelihoods, AIC/BIC, an AIC-based best-of-three selector, and normal QQ-plot residuals                                                                           |
| `probability`      | Bayes' theorem, joint/total probability, combinatorics (factorial, combinations, permutations), binomial and Poisson PMFs, and information theory (Shannon entropy, KL divergence, mutual information)                                                                                          |
| `copula`           | A bivariate Gaussian copula sampler with its own from-scratch normal CDF/inverse-CDF approximations, plus two rank-based dependence measures — Spearman's rho and Kendall's tau-b                                                                                                               |
| `correlation`      | Cholesky factorization of a correlation matrix for generating correlated normals, a separate (coarser) inverse-normal-CDF approximation (`probit`), and a distribution-quantile dispatcher (`icdf`) for QMC/Latin-hypercube sampling — no correlation-measurement function lives here           |

## `distributions`: twenty-six samplers behind one dispatcher, an unseeded RNG, and a gap at code 6

Every `sample_*` function draws directly from `std.random`'s `random_float64`/`randn_float64` with no seed argument anywhere in the module — unlike `discrete_event.queue_sim` elsewhere in this corpus, which takes an explicit seed and reproduces bit-for-bit on request, calling any function here twice with identical arguments gives two different draws. That makes a literal "here is the output" example misleading for most of these functions; what follows instead is what each family guarantees, plus the handful of genuinely deterministic edge cases.

| Family                                  | Functions                                                                                                                                                                                                                                                                 | Method                                                                                                                                                                                                                             |
| --------------------------------------- | ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- |
| Location-scale, closed-form inverse CDF | `sample_normal`, `sample_uniform`, `sample_lognormal`, `sample_exponential`, `sample_weibull`, `sample_laplace`, `sample_logistic`, `sample_gumbel`, `sample_rayleigh`, `sample_cauchy`, `sample_pareto`, `sample_arcsine`, `sample_halfnormal`                           | Transform a uniform (or normal) draw directly through the distribution's inverse CDF                                                                                                                                               |
| Shape-parameterized, rejection sampling | `sample_gamma_mt` (Marsaglia-Tsang), `sample_beta` (ratio of two gammas), `sample_student_t` (normal over sqrt of a scaled chi-squared draw), `sample_vonmises` (Best-Fisher rejection loop)                                                                              | An unbounded `while True` accept/reject loop — almost surely terminates, but has no iteration cap                                                                                                                                  |
| Bounded / project-specific              | `sample_triangular`, `sample_pert` (Beta-PERT via `sample_beta`)                                                                                                                                                                                                          | Closed-form inverse CDF (triangular) or a beta reparameterization (PERT)                                                                                                                                                           |
| Discrete counting                       | `sample_bernoulli`, `sample_poisson`, `sample_neg_binomial` (gamma-Poisson mixture), `sample_hypergeometric`                                                                                                                                                              | Direct simulation (hypergeometric draws without replacement in a loop) or inversion; `sample_poisson` switches to a normal approximation above `lam = 30`                                                                          |
| Universal dispatch                      | `sample_one(dt, p0, p1, p2, p3, dv, dp)`                                                                                                                                                                                                                                  | One `Int` type code selects among all of the above, plus `dt=4` (fixed constant, returns `p0`), `dt=13` (custom discrete via a cumulative-probability list `dp` and value list `dv`), and `dt=26` (categorical, returns the index) |
| SIMD4 batch path                        | `sample_normal_simd4`, `sample_uniform_simd4`, `sample_triangular_simd4`, `sample_lognormal_simd4`, `sample_exponential_simd4`, `sample_bernoulli_simd4`, `sample_gamma_simd4`, `sample_beta_simd4`, `sample_weibull_simd4`, `sample_pert_simd4`, plus `sample_one_simd4` | The same ten distributions, using an explicit `xorshift64` state threaded through the call (`mut s: UInt64`) instead of the stdlib RNG — four draws per call                                                                       |

**Two deterministic edge cases survive the randomness.** `sample_uniform(5.0, 5.0)` always returns exactly `5.0` — the formula is `lo + (hi - lo) * u`, and when `hi == lo` the `(hi - lo)` factor zeroes out whatever `u` was, independent of any RNG state. `sample_triangular(5.0, 5.0, 5.0)` and `sample_pert(5.0, 5.0, 5.0, 4.0)` are the same: both check `hi - lo == 0.0` and return `lo` directly before touching the RNG at all. These are the only calls in the module whose output does not depend on the draw.

**`dt = 6` is a genuine gap, not a hidden distribution.** `mojo/engine/types.mojo`'s own canonical `dist_type` table — the comment block that defines which integer means which distribution elsewhere in the engine — lists `5=beta` then jumps straight to `7=exponential`; `6` is not assigned to anything. `sample_one`'s dispatch chain matches that gap exactly: there is no `elif dt == 6` arm, so passing `dt=6` falls through to the same `return p0` fallback used for any unrecognized code (`dt=99` would behave identically) — a caller cannot tell "type 6" apart from a typo by the return value alone. Both the type table and the dispatcher agree with each other, so this is a deliberate numbering gap carried consistently through the codebase, not a mismatch between two disagreeing sources — but it does mean the module's own docstring, "all 27 native Mojo distributions," is counting slots `0`-`26` (27 of them), not 27 *distinct behaviors*: 26 codes are wired to a real sampler or discrete-value lookup, and one (`6`) is not.

## `distribution_fit`: MLE for two families, method-of-moments for a third, and an AIC upset

| Function                                                                            | What it does                                                                                                                                                                                                                                                                                            |
| ----------------------------------------------------------------------------------- | ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- |
| `fit_normal(data)`                                                                  | MLE: `[mean(data), population_std(data)]`. Empty input returns `[]` rather than dividing by zero.                                                                                                                                                                                                       |
| `fit_exponential(data)`                                                             | MLE: `rate = 1 / mean(data)`. Returns `1.0` if the mean is non-positive, `0.0` on empty input.                                                                                                                                                                                                          |
| `fit_uniform(data)`                                                                 | `[min(data), max(data)]` — the MLE support bounds for a uniform distribution.                                                                                                                                                                                                                           |
| `fit_weibull(data)`                                                                 | Method-of-moments approximation: coefficient of variation → shape (`1.2 / cv`, floored at `0.1`), then scale = `mean / Γ(1 + 1/shape)` using a quadratic gamma-function stand-in, `_gamma_approx`.                                                                                                      |
| `log_likelihood_normal(data, mean, std)` / `log_likelihood_exponential(data, rate)` | Closed-form log-likelihood for each family, evaluated at the given parameters.                                                                                                                                                                                                                          |
| `aic(log_likelihood, k)`                                                            | `2k - 2·LL`.                                                                                                                                                                                                                                                                                            |
| `bic(log_likelihood, k, n)`                                                         | `k·ln(n) - 2·LL`.                                                                                                                                                                                                                                                                                       |
| `best_fit(data)`                                                                    | Fits normal, exponential, and a uniform-with-support-check log-likelihood, computes AIC for each (`k=2, 1, 2` parameters respectively), and returns the index (`0`/`1`/`2`) of the lowest. Returns `-1` on empty input, guarding the `params[0]`/`params[1]` reads the three fits would otherwise need. |
| `qq_residuals(data, mean, std)`                                                     | Sorts `data`, then returns `sorted[i] - (mean + std·Φ⁻¹((i+0.5)/n))` per point, using its own inverse-erf-based `_normal_quantile`.                                                                                                                                                                     |

**Worked comparison — `best_fit` genuinely can pick the "wrong-looking" family.** For `data = [1.0, 2.0, 3.0, 4.0, 5.0]`, `fit_normal` gives `mean = 3.0`, `std = 1.4142135623730951` (population standard deviation, dividing by `n` not `n-1`), and that plugged into `log_likelihood_normal` gives `LL = -8.827560617423226`, so `aic(LL, 2) = 21.65512123484645`. `fit_exponential` gives `rate = 0.3333333333333333`, `LL = -10.49306144334055`, `aic(LL, 1) = 22.9861228866811`. `fit_uniform` gives `[1.0, 5.0]`; its own internal `_log_likelihood_uniform` (not exposed on its own) gives `LL = -6.931471805599453` (`= -n·ln(4)`), so `aic(LL, 2) = 17.862943611198908`. Uniform's AIC is the lowest of the three, so `best_fit([1.0, 2.0, 3.0, 4.0, 5.0])` returns `2`, not `0` — on five evenly-spaced points, a flat distribution over `[1, 5]` really is a better AIC fit than a bell curve centered at `3`, because AIC only rewards likelihood at the cost of parameter count, and evenly-spaced data has no tails for a normal distribution to explain away. Every number above was independently re-derived from the formulas in `mojo/engine/distribution_fit.mojo` in plain Python (see [Verification](#verification)).

**`fit_weibull`'s docstring says "iterative estimation"; the code takes one direct pass.** For the same five points, `cv = std/mean = 0.4714...`, `shape = 1.2/cv = 2.5455844122715705`, and `_gamma_approx(1 + 1/shape) = 0.8914274653601822` (against SciPy's/`math.gamma`'s true `Γ(1.3928371...) = 0.8876774286375584` — about a 0.42% relative error from the quadratic interpolation `_gamma_approx` uses for its `[1, 2)` range), giving `scale = mean / gamma_approx = 3.365388791097935`. There is no loop refining `shape` or `scale` against the data a second time; "iterative" in the docstring is closer to describing `_gamma_approx`'s own recursive range-reduction (it recurses to shift its argument into `[1, 2)` before applying the quadratic) than to the fit itself, which is a single closed-form CV-matching step.

## `probability`: Bayes, combinatorics, and information theory

| Function                                 | What it does                                                                                                                                                             |
| ---------------------------------------- | ------------------------------------------------------------------------------------------------------------------------------------------------------------------------ |
| `bayes(prior, likelihood, evidence)`     | `P(A\|B) = likelihood · prior / evidence`; returns `0.0` if `evidence` is within `1e-15` of zero.                                                                        |
| `joint_probability(p_a, p_b_given_a)`    | `P(A∩B) = p_b_given_a · p_a`.                                                                                                                                            |
| `total_probability(priors, likelihoods)` | `Σ priors[i]·likelihoods[i]` over `min(len(priors), len(likelihoods))` terms.                                                                                            |
| `factorial(n)`                           | Iterative `n!`; `1.0` for `n <= 1`.                                                                                                                                      |
| `combinations(n, k)`                     | `C(n,k)` via the multiplicative formula (using the smaller of `k`, `n-k`); `0.0` outside `0 <= k <= n`.                                                                  |
| `permutations(n, k)`                     | `P(n,k) = n!/(n-k)!`; `0.0` outside `0 <= k <= n`.                                                                                                                       |
| `binomial_pmf(k, n, p)`                  | `C(n,k)·pᵏ·(1-p)ⁿ⁻ᵏ`; `0.0` outside `0 <= k <= n`.                                                                                                                       |
| `poisson_pmf(k, lam)`                    | `λᵏ·e⁻λ/k!`; `0.0` for `k < 0`.                                                                                                                                          |
| `entropy(probs)`                         | Shannon entropy in **bits**: `-Σ p·log₂(p)`, skipping entries with `p <= 1e-15` (the `0·log 0 = 0` convention).                                                          |
| `kl_divergence(p, q)`                    | `Σ p·log₂(p/q)` in bits over the shared prefix of both lists; adds a `1e10` sentinel term per entry where `p > 0` but `q ≈ 0` (undefined KL) instead of producing `inf`. |
| `mutual_information(joint)`              | `I(X;Y) = ΣΣ p(x,y)·log₂(p(x,y)/(p(x)p(y)))` in bits, computing the marginals from the 2-D `joint` table itself.                                                         |

**The disease-test example in the opening code block is the classic illustration of why `bayes` matters**: even at 99% sensitivity, a 1% base rate and a 5% false-positive rate combine to an evidence term of `0.0594`, and `bayes(0.01, 0.99, 0.0594) = 0.16666...` — a positive result moves the odds from 1-in-100 to only 1-in-6, not to near-certainty. `total_probability([0.3, 0.7], [0.9, 0.2]) = 0.3·0.9 + 0.7·0.2 = 0.41` is the same law applied without collapsing to a single posterior.

**Entropy, KL, and mutual information all check out against their textbook extremes.** `entropy([0.5, 0.5]) = 1.0` bit exactly (a fair coin carries one bit of surprise) and `entropy([0.25, 0.25, 0.25, 0.25]) = 2.0` bits exactly (four equally likely outcomes, `log₂4`). `kl_divergence([0.5, 0.5], [0.9, 0.1]) = 0.7369655941662062` bits — the cost, in bits, of coding a fair-coin source with a model that expects it to land heads 90% of the time. For a joint table with independent margins (`[[0.25, 0.25], [0.25, 0.25]]`), `mutual_information` returns exactly `0.0`; for a dependent one concentrated on the diagonal (`[[0.4, 0.1], [0.1, 0.4]]`), it returns `0.2780719051126378` bits of shared information between the two variables. Every one of these was independently recomputed in Python from the formulas above and matches to full double precision.

## `copula` and `correlation`: dependence beyond a single number, in two unrelated files

**`copula.gaussian_copula_sample(n, rho)`** draws `n` bivariate pairs `(u1, u2)`, each marginally uniform on `[0, 1]`, whose *underlying* Gaussian latent variables are correlated at `rho`. It does this entirely by hand: its own `_box_muller_pair` (Box-Muller via a from-scratch Taylor-series `_sin_approx`/`_cos_approx`, not `std.math`'s), a hardcoded 2×2 Cholesky (`x2 = rho·z1 + sqrt(1-rho²)·z2`) rather than a call to `correlation.cholesky_lower`, and its own `norm_cdf_approx` to map the correlated normals back to `[0, 1]`. Because the RNG is unseeded, no single call's output is reproducible, but the *statistical* property is checkable: running the exact algorithm 200,000 times at `rho = 0.7` (Python's own RNG standing in for Mojo's, since the shape of the construction — not the specific stream — is what's being checked) gives a Spearman correlation of the resulting `u1, u2` pairs of `0.6827`, against the closed-form Gaussian-copula identity `6/π·arcsin(ρ/2) = 0.68291`, and a Kendall's tau of `0.4938` against `2/π·arcsin(ρ) = 0.49363` — both within the noise expected at that sample size, confirming the implementation reproduces the standard Gaussian copula rather than some other dependence structure.

**`copula.rank_correlation` and `copula.kendall_tau`** are ordinary Spearman and Kendall's tau-b, computed deterministically from data (no RNG involved) and independently verified against SciPy. For `x = [1, 2, 3, 4, 5]`, `y = [2, 1, 4, 3, 5]`: `rank_correlation(x, y) = 0.8` and `kendall_tau(x, y) = 0.6`, matching `scipy.stats.spearmanr`/`kendalltau` exactly. Both guard the same failure mode the same way — `n < 2` or mismatched lengths return `0.0` rather than reading past the shorter list's end (an uncatchable abort otherwise), per the source's own comments. Tied values are handled correctly too: `_rank` assigns the average rank across a tie block before either function runs Pearson-on-ranks (Spearman) or counts concordant/discordant pairs with a separate tie tally (Kendall's tau-b denominator), verified against SciPy on a tied input as well.

**`copula.norm_cdf_approx`/`norm_ppf_approx` and `correlation.probit` are two independent implementations of the same two functions, at two different precisions.** `norm_cdf_approx` is the Abramowitz & Stegun 26.2.17 approximation (stated max error `7.5e-8`); `norm_ppf_approx` is Peter Acklam's rational approximation, accurate to about `1e-9` — reproducing it in Python and comparing to `scipy.stats.norm.ppf` confirms errors of `1.6e-9` to `3.5e-10` across `p = 0.9, 0.975, 0.99`. `correlation.probit`, by contrast, uses the coarser Abramowitz & Stegun 26.2.23 rational approximation (stated max error `4.5e-4`) — the same comparison shows errors of `1.8e-4` to `4.4e-4` at the same `p` values, three to four orders of magnitude worse than `norm_ppf_approx`. Neither module imports the other's version; a caller who needs the more accurate inverse CDF and reaches for `correlation.probit` instead of `copula.norm_ppf_approx` gets a real, measurable accuracy loss, not a rounding difference.

**`correlation.cholesky_lower` and `correlation.icdf`** are the two building blocks for scenario generation. `cholesky_lower(corr, n)` factors an `n×n` correlation matrix into a lower triangular `L` with `L·Lᵀ = corr`; for `corr = [[1.0, 0.5, 0.3], [0.5, 1.0, 0.2], [0.3, 0.2, 1.0]]` it returns `[[1.0, 0, 0], [0.5, 0.8660254037844386, 0], [0.3, 0.0577350269189626, 0.9521904571390466]]`, matching NumPy's `linalg.cholesky` to full double precision. A non-positive-semidefinite `corr` does not crash it: a negative diagonal term is clamped to `0.0` before the square root (`sqrt(max(diag, 0.0))`) rather than producing a `NaN`. `sample_correlated_normal(L, n_vars)` then multiplies `L` by `n_vars` fresh `randn_float64()` draws to produce one vector of correlated standard normals — non-deterministic per call, like everything in `distributions`. `icdf(u, dt, p0, p1, p2)` is a separate, closed-form quantile dispatcher for QMC/Latin-hypercube sampling, covering the subset of `distributions.sample_one`'s type codes that have a closed-form inverse CDF — normal (via `probit`), uniform, triangular, lognormal, exponential, and Weibull, plus `dt=4` for a fixed constant — and deliberately omits the ones that don't, such as beta, gamma, and Poisson, which `sample_one` only reaches via rejection sampling or a series expansion.

{% hint style="warning" %}
**`cholesky_lower` has no guard against a `corr` shorter than `n`; `sample_correlated_normal` does.** Reading `cholesky_lower`'s body directly: it indexes `corr[i][i]` and `corr[i][j]` for every `i, j < n` with no bounds check at all — if the caller passes a `corr` with fewer than `n` rows, that is an unguarded, uncatchable out-of-bounds read that aborts the worker. `sample_correlated_normal`, by contrast, explicitly checks `len(L) < n_vars` and every row's length before touching `L[i][j]`, returning an empty list rather than aborting, and its own comment says why: "an unguarded `L[i][j]` read aborts the worker (uncatchable)." The runtime contract's `count_relations` entry for `cholesky_lower` (binding `corr`'s length to `n`) is exactly the enforcement point that has to catch a mismatch before the call reaches this kernel, since the kernel itself will not.
{% endhint %}

## Verification

All five `mojo/tests/libs/test_{distributions,distribution_fit,probability,copula,correlation}.mojo` files are the same auto-generated coverage-smoke pattern used across this corpus: every function is called once on a shared `0.5`-valued fixture (or `5`/`0.5` for the `Int`-typed distribution-type codes) and the test only checks that the result is finite (or, for list-returning functions, merely non-empty) — none of the five pin a specific expected number, and none of them exercise `cholesky_lower` or `sample_correlated_normal` at all (only `probit` and `icdf` are covered in `test_correlation.mojo`, since the auto-generator's shared fixture has no ready-made `List[List[Float64]]` value to hand `cholesky_lower`'s `corr` parameter).

Every worked number in this page — the Bayes disease-test posterior, the `best_fit` AIC comparison and its normal/exponential/uniform log-likelihoods, the `fit_weibull` shape/scale/gamma-approximation figures, the entropy/KL/mutual-information triple, the Spearman/Kendall pair on `[1,2,3,4,5]` vs. `[2,1,4,3,5]`, the Gaussian-copula rank-correlation identities, the `norm_ppf_approx`-vs-`probit` precision comparison, and the 3×3 Cholesky factorization — was independently re-derived from the real formulas in `mojo/engine/{distributions,distribution_fit,probability,copula,correlation}.mojo` using a standalone Python script (plain `math`, cross-checked against `scipy.stats` and `numpy.linalg` where a reference implementation exists), not lifted from a test file or asserted from memory. Values that depend on `exp`/`log`/`sqrt` may differ from a live Mojo run in the last one to a few representable digits — Mojo's own transcendental functions round a fraction of an ULP differently from CPython's, the same caveat noted elsewhere in this corpus — but no worked number here depends on that last-digit agreement to make its point.

**Next:** [Statistical inference](/time-series-and-statistics/statistical-inference.md) covers the question that comes after a distribution is fit — is an observed difference real, and how much data does it take to know — and [Discrete-event simulation](/optimization-and-simulation/discrete-event-simulation.md) is the one place elsewhere in this corpus where a random draw *is* seeded and reproducible, the opposite of every sampler in `distributions` here.


---

# Agent Instructions
This documentation is published with GitBook. GitBook is the documentation platform designed so that both humans and AI agents can read, navigate, and reason over technical content effectively. Learn more at gitbook.com.

## Querying This Documentation
If you need additional information that is not directly available in this page, you can query the documentation by asking a question.

Perform an HTTP GET request on the following URL with the `ask` and `goal` query parameters:

```
GET https://docs.algenta.ai/numerical-and-mathematical-computing/probability-distributions-and-dependence-kernels.md?ask=<question>&goal=<user_goal>
```

`ask` is the immediate question: it should be specific, self-contained, and written in natural language.
`goal` is what the user is ultimately trying to achieve, the reason they need the answer. Sharing it helps GitBook give you a better, more relevant answer. A goal is most helpful when it describes the outcome the user wants rather than restating the question. For example, with `ask=how do I create an API token`, a goal like `automate deployments from our CI pipeline` lets GitBook tailor the answer to that use case.

The response will contain a direct answer to the question and relevant excerpts and sources from the documentation.

Use this mechanism when the answer is not explicitly present in the current page, you need clarification or additional context, or you want to retrieve related documentation sections.