Use scipy.stats functions for defining distributions

[ryants]

Currently, priors/posteriors are calculated with our own code. We could use scipy and have numerous other distributions available, while not having to write/check the code ourselves

[wshlavacek]

Current state

Distribution handling is still entirely custom code:

• Sampling (pset.py): Uses raw np.random.normal / np.random.uniform with manual log-space transformations in sample_value().
• Prior calculation (algorithms.py, ln_prior()): Hand-rolled log-PDF formulas for normal and uniform priors. The normal prior omits the normalization constant (fine for MCMC acceptance ratios, but fragile if the code is reused elsewhere).

Unit tests added

Committed 29 regression tests in tests/test_distributions.py (393eb23) that lock in the current behavior of all 4 distribution types and all prior calculations, including verification against scipy.stats reference values. These will guard against regressions during the migration.

Implementation plan

Step 1: Replace sampling with scipy.stats distribution objects

In FreeParameter.__init__(), replace:

self._distribution = np.random.normal  # or np.random.uniform

with scipy.stats distribution objects that encapsulate both sampling and PDF evaluation:

• normal_var(mu, sigma) → scipy.stats.norm(loc=mu, scale=sigma)
• lognormal_var(mu, sigma) → sample via scipy.stats.norm(loc=mu, scale=sigma) in log10-space, then transform (scipy's lognorm uses natural log and has a non-obvious parameterization, so it's cleaner to keep the log10 transform explicit)
• uniform_var(lo, hi) → scipy.stats.uniform(loc=lo, scale=hi-lo)
• loguniform_var(lo, hi) → scipy.stats.uniform(loc=log10(lo), scale=log10(hi)-log10(lo)) in log10-space, then transform

Then sample_value() becomes self._distribution.rvs() with the appropriate transform.

Step 2: Replace ln_prior() with distribution logpdf

Store the scipy.stats distribution on each FreeParameter and use dist.logpdf() in BayesianAlgorithm.ln_prior(). This eliminates the hand-rolled formulas and the separate prior dict of tuples.

Note: the current normal prior omits the normalization constant -0.5*log(2*pi*sigma^2). Switching to scipy.stats.norm.logpdf() will include it, which changes absolute values but not MCMC acceptance ratios (the constant cancels in differences). The regression tests verify this via the TestPriorDifferenceForMCMC test class.

Step 3: Add new prior distribution types

Once distributions are scipy.stats objects, adding new prior types becomes trivial — just map a config keyword to a scipy.stats distribution. Candidates:

• Gamma (scipy.stats.gamma) — common prior for rate parameters
• Beta (scipy.stats.beta) — natural prior for probabilities / fractions
• Cauchy / half-Cauchy (scipy.stats.cauchy, scipy.stats.halfcauchy) — heavy-tailed priors, popular in hierarchical Bayesian models
• Exponential (scipy.stats.expon) — memoryless prior for scale parameters
• Truncated normal (scipy.stats.truncnorm) — bounded normal, useful when parameters have physical constraints

Each new distribution type requires only: a new config keyword, a line mapping it to a scipy.stats object, and a parser entry.

Benefits

• Expanded Bayesian capabilities: Users can specify informative priors from a much richer family of distributions, improving the quality of posterior inference for domain-specific problems.
• Better maintainability: Replaces ~50 lines of hand-rolled distribution math with well-tested, numerically stable scipy.stats calls. Fewer places for bugs to hide.
• Consistency: Sampling and PDF evaluation use the same distribution object, eliminating the risk of parameterization mismatches between the two code paths.
• Numerical stability: scipy.stats.logpdf() handles edge cases (extreme values, underflow) more robustly than the current manual formulas.