mathematicskit.probability#

mathematicskit.probability hub · History · Examples · API reference

Discrete and continuous distributions (binomial, Poisson, geometric, uniform, exponential, normal, gamma) built on scipy.stats, with mathematicskit’s own moment generating functions; Monte Carlo integration with variance reduction (importance sampling, control variates); Law of Large Numbers and Central Limit Theorem simulation; and discrete-time Markov chains (stationary distributions, absorption probabilities).

mathematicskit.probability: probability theory, built directly on scipy.stats.

Discrete and continuous distribution classes (binomial, Poisson, geometric, uniform, exponential, normal, gamma, beta) wrapping scipy.stats, with mathematicskit’s own moment generating function added on top; law-of-large-numbers and central-limit-theorem simulation via numpy.random sampling; Monte Carlo integration with variance reduction (importance sampling, control variates, hand-rolled – no direct scipy equivalent); and discrete-time Markov chains (stationary distribution via numpy.linalg.eig/power iteration, absorption probabilities and expected absorption time via numpy.linalg.solve). Also: Buffon’s needle, the St. Petersburg game, Beta-binomial Bayesian updating, Chebyshev’s inequality, Galton-Watson branching processes, random walks and Brownian motion, Erlang’s loss formula, continuous-time Markov chains via scipy.linalg.expm(), Markov chain Monte Carlo (Metropolis-Hastings and the Gibbs sampler), and stochastic differential equations (Euler-Maruyama, geometric Brownian motion, Ornstein-Uhlenbeck).

class mathematicskit.probability.Beta(alpha, beta)[source]#

Bases: ContinuousDistribution

Beta distribution with shape parameters \(\alpha\) and \(\beta\) on \([0, 1]\).

\(f(x) = \dfrac{x^{\alpha-1}(1-x)^{\beta-1}}{B(\alpha, \beta)}\), via scipy.stats.beta. MGF: \(M(t) = {}_1F_1(\alpha; \alpha+\beta; t)\), Kummer’s confluent hypergeometric function (scipy.special.hyp1f1()). It is the conjugate prior of the binomial, which makes it the natural home of Bayes’s 1763 problem (see mathematicskit.probability.systems.bayes). See DeGroot & Schervish, Probability and Statistics, 4th ed., Sec. 5.8.

Parameters:

  • alpha (float) – Shape parameters, both > 0.

  • beta (float) – Shape parameters, both > 0.

Examples

>>> b = Beta(alpha=2.0, beta=3.0)
>>> round(b.mean, 4)
0.4
>>> round(b.variance, 4)
0.04
mgf(t)[source]#

Moment generating function, E[e^{tX}] (closed form; not exposed by scipy.stats).

pdf(x)[source]#

Probability density function.

class mathematicskit.probability.Binomial(n, p)[source]#

Bases: DiscreteDistribution

Binomial distribution: number of successes in n i.i.d. Bernoulli(p) trials.

\(P(X=k) = \binom{n}{k}p^k(1-p)^{n-k}\), via scipy.stats.binom. MGF: \(M(t) = (1-p+pe^t)^n\). See DeGroot & Schervish, Probability and Statistics, 4th ed., Sec. 5.4.

Parameters:

  • n (int) – Number of trials.

  • p (float) – Success probability, 0 <= p <= 1.

Examples

>>> b = Binomial(n=10, p=0.3)
>>> round(b.mean, 4)
3.0
>>> round(b.variance, 4)
2.1
>>> round(float(b.mgf(0.0)), 10)
1.0
mgf(t)[source]#

Moment generating function, E[e^{tX}] (closed form; not exposed by scipy.stats).

pmf(k)[source]#

Probability mass function, P(X = k).

class mathematicskit.probability.BranchingProcessResult(generation_sizes, extinct_fraction)[source]#

Bases: object

Container for simulated Galton-Watson generation sizes.

Parameters:
extinct_fraction: float#

Fraction of runs whose population reached zero within the simulated generations.

Type:

float

generation_sizes: ndarray#

Population of each run at each generation (column 0 is the founders).

Type:

ndarray, shape (n_runs, n_generations + 1)

class mathematicskit.probability.BrownianMotionResult(times, paths)[source]#

Bases: object

Container for sampled Brownian-motion paths.

Parameters:
paths: ndarray#

Sampled paths W(t), each starting at W(0) = 0.

Type:

ndarray, shape (n_paths, n_steps + 1)

times: ndarray#

Time grid, starting at 0.

Type:

ndarray, shape (n_steps + 1,)

class mathematicskit.probability.BuffonNeedleResult(pi_estimate, crossing_fraction, exact_crossing_probability, n_drops)[source]#

Bases: object

Container for a Buffon’s-needle simulation.

Parameters:

  • pi_estimate (float)

  • crossing_fraction (float)

  • exact_crossing_probability (float)

  • n_drops (int)

crossing_fraction: float#

Fraction of dropped needles that crossed a line.

Type:

float

exact_crossing_probability: float#

The exact crossing probability 2 * length / (pi * spacing) for a short needle.

Type:

float

n_drops: int#

Number of needles dropped.

Type:

int

pi_estimate: float#

The estimate of pi implied by the observed crossing fraction.

Type:

float

class mathematicskit.probability.ContinuousDistribution[source]#

Bases: ABC

Common base for scipy.stats-backed continuous distribution wrappers.

cdf(x)[source]#

Cumulative distribution function, P(X <= x).

property mean: float#

E[X].

Type:

float

abstractmethod mgf(t)[source]#

Moment generating function, E[e^{tX}] (closed form; not exposed by scipy.stats).

abstractmethod pdf(x)[source]#

Probability density function.

ppf(q)[source]#

Quantile function (inverse CDF), P(X <= ppf(q)) = q.

sample(size=1, seed=None)[source]#

Draw samples via scipy.stats’ rvs.

Parameters:
Returns:

ndarray

property std: float#

sqrt(Var(X)).

Type:

float

property variance: float#

Var(X).

Type:

float

class mathematicskit.probability.DiscreteDistribution[source]#

Bases: ABC

Common base for scipy.stats-backed discrete distribution wrappers.

cdf(k)[source]#

Cumulative distribution function, P(X <= k).

property mean: float#

E[X].

Type:

float

abstractmethod mgf(t)[source]#

Moment generating function, E[e^{tX}] (closed form; not exposed by scipy.stats).

abstractmethod pmf(k)[source]#

Probability mass function, P(X = k).

ppf(q)[source]#

Quantile function (inverse CDF), P(X <= ppf(q)) = q.

sample(size=1, seed=None)[source]#

Draw samples via scipy.stats’ rvs.

Parameters:
Returns:

ndarray

property std: float#

sqrt(Var(X)).

Type:

float

property variance: float#

Var(X).

Type:

float

class mathematicskit.probability.Exponential(rate)[source]#

Bases: ContinuousDistribution

Exponential distribution: waiting time between Poisson events of rate rate.

\(f(x) = \text{rate}\,e^{-\text{rate}\,x}\), \(x \geq 0\), via scipy.stats.expon. MGF: \(M(t) = \dfrac{\text{rate}}{\text{rate}-t}\) for \(t < \text{rate}\). See DeGroot & Schervish, Probability and Statistics, 4th ed., Sec. 5.7.

Parameters:

rate (float) – Rate parameter \(\lambda > 0\).

Examples

>>> e = Exponential(rate=2.0)
>>> round(e.mean, 4)
0.5
>>> round(e.variance, 4)
0.25
mgf(t)[source]#

Moment generating function, E[e^{tX}] (closed form; not exposed by scipy.stats).

pdf(x)[source]#

Probability density function.

class mathematicskit.probability.Gamma(shape, rate)[source]#

Bases: ContinuousDistribution

Gamma distribution with shape k and rate \(\beta\).

\(f(x) = \dfrac{\beta^k}{\Gamma(k)}x^{k-1}e^{-\beta x}\), \(x \geq 0\), via scipy.stats.gamma. MGF: \(M(t) = \left(\dfrac{\beta}{\beta-t}\right)^k\) for \(t < \beta\). See DeGroot & Schervish, Probability and Statistics, 4th ed., Sec. 5.8.

Parameters:

  • shape (float) – Shape parameter k > 0.

  • rate (float) – Rate parameter \(\beta > 0\).

Examples

>>> g = Gamma(shape=2.0, rate=1.0)
>>> round(g.mean, 4)
2.0
>>> round(g.variance, 4)
2.0
mgf(t)[source]#

Moment generating function, E[e^{tX}] (closed form; not exposed by scipy.stats).

pdf(x)[source]#

Probability density function.

class mathematicskit.probability.Geometric(p)[source]#

Bases: DiscreteDistribution

Geometric distribution: number of trials up to and including the first success.

\(P(X=k) = (1-p)^{k-1}p\), k = 1, 2, \dots (scipy’s “number of trials” convention), via scipy.stats.geom. MGF: \(M(t) = \dfrac{pe^t}{1-(1-p)e^t}\), for \(t < -\ln(1-p)\). See DeGroot & Schervish, Probability and Statistics, 4th ed., Sec. 5.3.

Parameters:

p (float) – Success probability per trial, 0 < p <= 1.

Examples

>>> g = Geometric(p=0.25)
>>> round(g.mean, 4)
4.0
mgf(t)[source]#

Moment generating function, E[e^{tX}] (closed form; not exposed by scipy.stats).

pmf(k)[source]#

Probability mass function, P(X = k).

class mathematicskit.probability.MCMCResult(samples, acceptance_rate, method='')[source]#

Bases: object

Container for a Markov chain Monte Carlo run.

Parameters:
acceptance_rate: float#

Fraction of proposals accepted (1.0 for a Gibbs sampler, which never rejects).

Type:

float

method: str = ''#

"metropolis_hastings" or "gibbs".

Type:

str

samples: ndarray#

The chain’s states after burn-in, one row per step.

Type:

ndarray, shape (n_samples, dim)

class mathematicskit.probability.MarkovChain(transition_matrix)[source]#

Bases: object

A discrete-time, finite-state Markov chain given its transition matrix.

Parameters:

transition_matrix (ndarray) – Row-stochastic matrix: transition_matrix[i, j] is \(P(X_{t+1}=j \mid X_t=i)\), so each row sums to 1.

Examples

>>> import numpy as np
>>> # A simple 2-state weather chain: sunny/rainy.
>>> P = np.array([[0.9, 0.1], [0.5, 0.5]])
>>> chain = MarkovChain(P)
>>> pi = chain.stationary_distribution()
>>> np.allclose(pi, [5.0 / 6.0, 1.0 / 6.0], atol=1e-8)
True
absorption_probabilities(transient, absorbing)[source]#

Probability of ending in each absorbing state, from each transient state.

Partitioning P (after reordering) into transient-to-transient block \(Q\) and transient-to-absorbing block \(R\), the absorption probabilities are \(B = (I - Q)^{-1} R\), solved via numpy.linalg.solve() (never an explicit matrix inverse). See Grinstead & Snell, Introduction to Probability, 2nd ed., Theorem 11.16.

Parameters:

  • transient (Sequence[int]) – Indices of the transient states.

  • absorbing (Sequence[int]) – Indices of the absorbing states (rows with a single 1 on the diagonal).

Return type:

ndarray

Returns:

ndarray, shape (len(transient), len(absorbing)) – B[i, j]: probability of eventual absorption in absorbing[j], starting from transient[i].

Examples

>>> import numpy as np
>>> # Gambler's ruin, capital 1..3, absorbing at 0 and 4, fair coin (p=0.5).
>>> P = np.array([
...     [1.0, 0.0, 0.0, 0.0, 0.0],
...     [0.5, 0.0, 0.5, 0.0, 0.0],
...     [0.0, 0.5, 0.0, 0.5, 0.0],
...     [0.0, 0.0, 0.5, 0.0, 0.5],
...     [0.0, 0.0, 0.0, 0.0, 1.0],
... ])
>>> chain = MarkovChain(P)
>>> B = chain.absorption_probabilities(transient=[1, 2, 3], absorbing=[0, 4])
>>> np.allclose(B[:, 1], [0.25, 0.5, 0.75], atol=1e-8)  # P(reach 4 | start at i) = i/4
True
expected_steps_to_absorption(transient)[source]#

Expected number of steps before absorption, from each transient state.

\(t = (I-Q)^{-1}\mathbf{1}\), via numpy.linalg.solve(). See Grinstead & Snell, Introduction to Probability, 2nd ed., Theorem 11.14.

Parameters:

transient (Sequence[int]) – Indices of the transient states.

Return type:

ndarray

Returns:

ndarray, shape (len(transient),)

Examples

>>> import numpy as np
>>> P = np.array([
...     [1.0, 0.0, 0.0, 0.0, 0.0],
...     [0.5, 0.0, 0.5, 0.0, 0.0],
...     [0.0, 0.5, 0.0, 0.5, 0.0],
...     [0.0, 0.0, 0.5, 0.0, 0.5],
...     [0.0, 0.0, 0.0, 0.0, 1.0],
... ])
>>> chain = MarkovChain(P)
>>> t = chain.expected_steps_to_absorption(transient=[1, 2, 3])
>>> np.allclose(t, [3.0, 4.0, 3.0], atol=1e-8)  # i * (4 - i) for i=1,2,3
True
stationary_distribution()[source]#

Stationary distribution \(\pi\) with \(\pi P = \pi\), \(\sum_i \pi_i = 1\).

\(\pi\) is the left eigenvector of P for eigenvalue 1, equivalently the right eigenvector of \(P^T\) – computed via numpy.linalg.eig(). Only well-defined (unique) for an irreducible chain; the eigenvalue closest to 1 is used, which for a reducible/periodic chain may not be exactly 1. See Grinstead & Snell, Introduction to Probability, 2nd ed., Theorem 11.6.

Return type:

ndarray

Returns:

ndarray, shape (n,)

stationary_distribution_power_iteration(n_iter=10000, tol=1e-14)[source]#

Stationary distribution via power iteration: \(\pi_{k+1} = \pi_k P\).

Starting from the uniform distribution, repeatedly applies the transition matrix; for an irreducible, aperiodic (“ergodic”) chain this converges to the same stationary \(\pi\) found by stationary_distribution(), illustrating the chain “forgetting” its start. See Grinstead & Snell, Introduction to Probability, 2nd ed., Theorem 11.7.

Parameters:

  • n_iter (int)

  • tol (float) – Convergence tolerance on the total-variation change per step.

Return type:

ndarray

Returns:

ndarray, shape (n,)

Examples

>>> import numpy as np
>>> P = np.array([[0.9, 0.1], [0.5, 0.5]])
>>> chain = MarkovChain(P)
>>> pi = chain.stationary_distribution_power_iteration()
>>> np.allclose(pi, [5.0 / 6.0, 1.0 / 6.0], atol=1e-6)
True
class mathematicskit.probability.MonteCarloResult(estimate, std_error, n_samples=0, method='')[source]#

Bases: object

Container for the output of a Monte Carlo integral estimate.

Parameters:
estimate: float#

The estimated integral/expectation.

Type:

float

method: str = ''#

e.g. "plain", "importance_sampling", "control_variates".

Type:

str

n_samples: int = 0#

Number of Monte Carlo samples used.

Type:

int

std_error: float#

Estimated standard error of estimate (std(samples) / sqrt(n), scaled appropriately by the method).

Type:

float

class mathematicskit.probability.Normal(mu, sigma)[source]#

Bases: ContinuousDistribution

Normal (Gaussian) distribution \(\mathcal{N}(\mu, \sigma^2)\).

Via scipy.stats.norm. MGF: \(M(t) = \exp(\mu t + \tfrac12\sigma^2 t^2)\). See DeGroot & Schervish, Probability and Statistics, 4th ed., Sec. 5.9.

Parameters:

  • mu (float) – Mean.

  • sigma (float) – Standard deviation, sigma > 0.

Examples

>>> n = Normal(mu=1.0, sigma=2.0)
>>> n.mean, n.variance
(1.0, 4.0)
mgf(t)[source]#

Moment generating function, E[e^{tX}] (closed form; not exposed by scipy.stats).

pdf(x)[source]#

Probability density function.

class mathematicskit.probability.Poisson(mu)[source]#

Bases: DiscreteDistribution

Poisson distribution: count of events in a fixed interval at rate mu.

\(P(X=k) = e^{-\mu}\mu^k/k!\), via scipy.stats.poisson. MGF: \(M(t) = e^{\mu(e^t - 1)}\). See DeGroot & Schervish, Probability and Statistics, 4th ed., Sec. 5.5.

Parameters:

mu (float) – Rate (mean number of events), mu > 0.

Examples

>>> p = Poisson(mu=4.0)
>>> p.mean == p.variance == 4.0
True
mgf(t)[source]#

Moment generating function, E[e^{tX}] (closed form; not exposed by scipy.stats).

pmf(k)[source]#

Probability mass function, P(X = k).

class mathematicskit.probability.SDEResult(times, paths, wiener, method='')[source]#

Bases: object

Container for sampled paths of a stochastic differential equation.

Parameters:
method: str = ''#

e.g. "euler_maruyama" or "exact".

Type:

str

paths: ndarray#

Sampled solution paths X(t).

Type:

ndarray, shape (n_paths, n_steps + 1)

times: ndarray#

Time grid, starting at 0.

Type:

ndarray, shape (n_steps + 1,)

wiener: ndarray#

The Brownian paths W(t) that drove them, so a numerical path can be compared with the exact solution driven by the same noise.

Type:

ndarray, shape (n_paths, n_steps + 1)

class mathematicskit.probability.TailBoundResult(k, bound, exact)[source]#

Bases: object

Container comparing Chebyshev’s bound with the exact two-sided tail probability.

Parameters:
bound: ndarray#

Chebyshev’s bound min(1, 1 / k**2).

Type:

ndarray

exact: ndarray#

The exact P(|X - mean| >= k * std) from the distribution’s CDF.

Type:

ndarray

k: ndarray#

Distances from the mean, in standard deviations.

Type:

ndarray

class mathematicskit.probability.Uniform(a, b)[source]#

Bases: ContinuousDistribution

Continuous uniform distribution on \([a, b]\).

\(f(x) = 1/(b-a)\) on \([a,b]\), via scipy.stats.uniform. MGF: \(M(t) = \dfrac{e^{tb}-e^{ta}}{t(b-a)}\) for \(t \neq 0\), \(M(0)=1\). See DeGroot & Schervish, Probability and Statistics, 4th ed., Sec. 5.6.

Parameters:

  • a (float) – Interval endpoints, a < b.

  • b (float) – Interval endpoints, a < b.

Examples

>>> u = Uniform(a=0.0, b=2.0)
>>> round(u.mean, 4)
1.0
>>> round(u.variance, 6)
0.333333
mgf(t)[source]#

Moment generating function, E[e^{tX}] (closed form; not exposed by scipy.stats).

pdf(x)[source]#

Probability density function.

mathematicskit.probability.beta_binomial_posterior(successes, trials, prior_alpha=1.0, prior_beta=1.0)[source]#

Posterior of a binomial success probability under a conjugate Beta prior.

By Bayes’s theorem, \(\pi(p \mid k) \propto p^k (1-p)^{n-k}\, \pi(p)\), and for a \(\mathrm{Beta}(\alpha, \beta)\) prior this is

\[p \mid k \sim \mathrm{Beta}(\alpha + k,\; \beta + n - k).\]
Parameters:

  • successes (int) – Observed successes \(k\).

  • trials (int) – Number of trials \(n \geq k\).

  • prior_alpha (float) – Prior shape parameters; the default (1, 1) is Bayes’s uniform prior.

  • prior_beta (float) – Prior shape parameters; the default (1, 1) is Bayes’s uniform prior.

Return type:

Beta

Returns:

Beta – The posterior distribution.

Examples

>>> post = beta_binomial_posterior(7, 10)
>>> post.alpha, post.beta
(8.0, 4.0)
>>> round(post.mean, 4)
0.6667
mathematicskit.probability.brownian_motion(n_paths=1, n_steps=1000, t_max=1.0, seed=0)[source]#

Sample standard Brownian motion on \([0, t_{\max}]\) on a uniform grid.

Each increment over a step \(\Delta t\) is an independent \(\mathcal N(0, \Delta t)\) draw, so \(W(t) \sim \mathcal N(0, t)\) exactly at every grid time.

Parameters:
Return type:

BrownianMotionResult

Returns:

BrownianMotionResult

Examples

>>> result = brownian_motion(n_paths=20000, n_steps=100, t_max=2.0, seed=0)
>>> result.paths.shape
(20000, 101)
>>> bool(abs(result.paths[:, -1].var() - 2.0) < 0.1)
True
mathematicskit.probability.buffon_needle(n_drops=100000, length=1.0, spacing=1.0, seed=0)[source]#

Drop n_drops needles on a ruled floor and estimate \(\pi\) from the crossings.

Each drop places the needle’s center a uniform distance \(y \in [0, d/2]\) from the nearest line and points it in a uniformly random direction. The needle crosses that line when \(y \leq \tfrac{\ell}{2}\sin\theta\), which happens with probability

\[P(\text{cross}) = \frac{2\ell}{\pi d}, \qquad \ell \leq d,\]

so \(\hat\pi = 2\ell n / (d \cdot \#\text{crossings})\). The random direction comes from normalizing a standard 2D Gaussian vector, so the simulation never uses \(\pi\) itself. See G.-L. Leclerc de Buffon, Essai d’arithmétique morale (1777), Sec. 23.

Parameters:

  • n_drops (int) – Number of needles dropped.

  • length (float) – Needle length \(\ell\).

  • spacing (float) – Distance \(d\) between the lines; must satisfy length <= spacing.

  • seed (int)

Return type:

BuffonNeedleResult

Returns:

BuffonNeedleResult

Examples

>>> result = buffon_needle(n_drops=400000, seed=0)
>>> bool(abs(result.pi_estimate - 3.14159) < 0.03)
True
>>> round(result.exact_crossing_probability, 4)
0.6366
mathematicskit.probability.central_limit_theorem_sample_means(dist, n, n_trials=2000, seed=0)[source]#

Standardized sample means from repeated trials, illustrating the CLT.

Draws n_trials independent samples of size n from dist, and returns the standardized sample mean \(Z = \dfrac{\bar X_n - E[X]}{\sigma/\sqrt n}\) for each trial; by the CLT, this converges in distribution to \(\mathcal N(0, 1)\) as \(n \to \infty\), regardless of dist’s own shape – verify with e.g. scipy.stats.kstest(result, "norm") (a two-sided goodness-of-fit test against the standard normal CDF). See DeGroot & Schervish, Probability and Statistics, 4th ed., Theorem 6.3.2.

Parameters:
Return type:

ndarray

Returns:

ndarray, shape (n_trials,)

Examples

>>> from mathematicskit.probability.systems.discrete import Poisson
>>> z = central_limit_theorem_sample_means(Poisson(mu=3.0), n=200, n_trials=5000, seed=0)
>>> abs(float(z.mean())) < 0.1
True
>>> abs(float(z.std()) - 1.0) < 0.1
True
mathematicskit.probability.chain_effective_sample_size(chain)[source]#

Effective number of independent draws in a correlated Markov chain, \(n/\tau\).

Unlike effective_sample_size() (for importance weights), this measures how much a chain’s autocorrelation costs: the mean of chain is as precise as the mean of this many independent draws.

Parameters:

chain (ndarray)

Return type:

float

Returns:

float

Examples

>>> import numpy as np
>>> chain = np.repeat(np.random.default_rng(1).standard_normal(2000), 10)  # each draw held for 10 steps
>>> bool(5 < 20000 / chain_effective_sample_size(chain) < 15)
True
mathematicskit.probability.chebyshev_tail(dist, k)[source]#

Chebyshev’s bound and the exact two-sided tail \(P(|X - \mu| \geq k\sigma)\).

\[P(|X - \mu| \geq k\sigma) \leq \frac{1}{k^2}.\]

The bound holds for every distribution with finite variance, which is both its strength and why it is loose for any particular one: for the normal distribution at \(k = 2\) it gives 0.25 against an exact 0.0455.

Parameters:
Return type:

TailBoundResult

Returns:

TailBoundResult

Examples

>>> from mathematicskit.probability.systems.continuous import Normal
>>> result = chebyshev_tail(Normal(mu=0.0, sigma=1.0), [2.0])
>>> float(result.bound[0]), round(float(result.exact[0]), 4)
(0.25, 0.0455)
mathematicskit.probability.control_variates_integrate(f, control, control_mean, a, b, n=10000, seed=0)[source]#

Control-variate Monte Carlo estimate of \(\int_a^b f(x)\,dx\).

Uses a correlated function control with known expectation control_mean under \(U(a,b)\) to reduce variance: estimates \((b-a)\,E[f(X) - c\,(g(X) - E[g(X)])]\), which has the same expectation as plain Monte Carlo but lower variance for the sample-covariance-optimal \(c = \mathrm{Cov}(f(X), g(X))/\mathrm{Var}(g(X))\). See Robert & Casella, Monte Carlo Statistical Methods, 2nd ed., Ch. 4.1.

Parameters:
Return type:

MonteCarloResult

Returns:

MonteCarloResult

Examples

>>> import numpy as np
>>> # Integrate exp(x) over [0, 1] = e - 1, using g(x) = x (E[X] = 0.5) as a control.
>>> result = control_variates_integrate(np.exp, lambda x: x, control_mean=0.5, a=0.0, b=1.0, n=200000, seed=0)
>>> abs(result.estimate - (np.e - 1.0)) < 0.01
True
mathematicskit.probability.ctmc_stationary_distribution(generator)[source]#

Stationary distribution \(\pi\) of an irreducible chain, solving \(\pi Q = 0\), \(\sum_i \pi_i = 1\).

Solved as a least-squares system with the normalization appended to \(Q^T\) (numpy.linalg.lstsq()); the solution is exact for an irreducible chain.

Parameters:

generator (array-like, shape (n, n))

Return type:

ndarray

Returns:

ndarray, shape (n,)

Examples

>>> ctmc_stationary_distribution([[-1.0, 1.0], [3.0, -3.0]]).round(4)
array([0.75, 0.25])
mathematicskit.probability.ctmc_transition_matrix(generator, t)[source]#

Transition matrix \(P(t) = e^{tQ}\) of a continuous-time Markov chain.

Parameters:

  • generator (array-like, shape (n, n)) – Generator matrix \(Q\).

  • t (float) – Elapsed time, t >= 0.

Return type:

ndarray

Returns:

ndarray, shape (n, n) – P[i, j] is the probability of being in state j at time t having started in state i.

Examples

>>> # Two states, rate 1 from 0 to 1 and rate 3 back: P_01(t) = (1 - e^{-4t}) / 4.
>>> P = ctmc_transition_matrix([[-1.0, 1.0], [3.0, -3.0]], t=0.5)
>>> bool(np.isclose(P[0, 1], (1 - np.exp(-2.0)) / 4))
True
mathematicskit.probability.effective_sample_size(weights)[source]#

Effective sample size of a set of (unnormalized) importance weights.

\(\mathrm{ESS} = \dfrac{\left(\sum_i w_i\right)^2}{\sum_i w_i^2}\), ranging from 1 (all weight on one sample – the proposal is a poor match for the target) to n (uniform weights – as good as plain Monte Carlo). A standard diagnostic for whether an importance sampler’s proposal distribution is well-matched to the integrand. See Robert & Casella, Monte Carlo Statistical Methods, 2nd ed., Ch. 3.3.2.

Parameters:

weights (ndarray) – Non-negative importance weights (e.g. f(x)/q(x) from importance_sampling_integrate()).

Return type:

float

Returns:

float

Examples

>>> import numpy as np
>>> round(effective_sample_size(np.ones(100)), 4)
100.0
>>> round(effective_sample_size(np.array([1.0, 0.0, 0.0, 0.0])), 4)
1.0
mathematicskit.probability.erlang_b(offered_load, servers)[source]#

Erlang B blocking probability for servers lines under offered_load erlangs.

Parameters:

  • offered_load (float) – Offered traffic \(A = \lambda / \mu \geq 0\), in erlangs.

  • servers (int) – Number of lines \(c \geq 0\).

Return type:

float

Returns:

float – The probability that an arriving call is blocked.

Examples

>>> round(erlang_b(1.0, 1), 4)
0.5
>>> round(erlang_b(2.0, 2), 4)
0.4
mathematicskit.probability.euler_maruyama(drift, diffusion, x0, t_max=1.0, n_steps=1000, n_paths=1, seed=0, dw=None)[source]#

Integrate the scalar Itô SDE \(dX = a(X, t)\,dt + b(X, t)\,dW\) by the Euler-Maruyama scheme.

Parameters:

  • drift (Callable[[ndarray, float], ndarray]) – drift(x, t) -> ndarray, vectorized over the array of paths x.

  • diffusion (Callable[[ndarray, float], ndarray]) – diffusion(x, t) -> ndarray, vectorized likewise.

  • x0 (float) – Initial value, shared by every path.

  • t_max (float)

  • n_steps (int)

  • n_paths (int)

  • seed (int) – Seeds the Brownian increments (ignored when dw is given).

  • dw (ndarray | None) – Brownian increments to use instead of fresh ones, for example a coarsened copy of a finer path’s increments, so that two step sizes can be compared on the same Brownian path.

Return type:

SDEResult

Returns:

SDEResult

Examples

>>> import numpy as np
>>> # dX = -X dt + dW from X(0) = 0 has Var X(t) = (1 - exp(-2t)) / 2.
>>> result = euler_maruyama(lambda x, t: -x, lambda x, t: np.ones_like(x), 0.0, t_max=1.0, n_steps=200, n_paths=20000)
>>> bool(abs(result.paths[:, -1].var() - (1 - np.exp(-2.0)) / 2) < 0.02)
True
mathematicskit.probability.galton_watson_extinction_probability(offspring_pmf, tol=1e-10, max_iter=1000)[source]#

Probability that a Galton-Watson family started by one individual eventually dies out.

Iterates \(q_{n+1} = G(q_n)\) from \(q_0 = 0\); each iterate is the probability of extinction by generation \(n\). In the subcritical and critical cases (\(m \leq 1\)) the answer is 1 by the Galton-Watson theorem and is returned directly, since the critical iteration converges only like \(1/n\).

Parameters:

  • offspring_pmf (array-like of float) – offspring_pmf[j] is the probability of exactly j children.

  • tol (float) – Stop when successive iterates differ by less than tol.

  • max_iter (int)

Return type:

float

Returns:

float – The extinction probability \(q\).

Examples

>>> # G(s) = 1/4 + s/4 + s^2/2, whose smallest fixed point is 1/2.
>>> round(galton_watson_extinction_probability([0.25, 0.25, 0.5]), 8)
0.5
mathematicskit.probability.galton_watson_simulate(offspring_pmf, n_generations, n_runs=1000, founders=1, seed=0)[source]#

Simulate independent Galton-Watson processes.

The children of the \(Z_n\) individuals in generation \(n\) are counted with one multinomial draw of size \(Z_n\) over the offspring classes, so \(Z_{n+1} = \sum_j j\,N_j\).

Parameters:

  • offspring_pmf (array-like of float)

  • n_generations (int)

  • n_runs (int) – Number of independent families.

  • founders (int) – Population of generation 0.

  • seed (int)

Return type:

BranchingProcessResult

Returns:

BranchingProcessResult

Examples

>>> result = galton_watson_simulate([0.25, 0.25, 0.5], n_generations=30, n_runs=4000, seed=0)
>>> bool(abs(result.extinct_fraction - 0.5) < 0.03)
True
mathematicskit.probability.geometric_brownian_motion(mu, sigma, x0=1.0, t_max=1.0, n_steps=1000, n_paths=1, seed=0, method='exact')[source]#

Sample geometric Brownian motion \(dX = \mu X\,dt + \sigma X\,dW\).

Itô’s formula applied to \(\log X\) gives the exact solution

\[X_t = X_0 \exp\!\left(\left(\mu - \tfrac12\sigma^2\right)t + \sigma W_t\right),\]

whose \(-\tfrac12\sigma^2\) is the Itô correction that ordinary calculus would miss. method="euler_maruyama" integrates the SDE instead; with the same seed both methods use the same Brownian path, so their difference is the scheme’s pathwise error.

Parameters:
Return type:

SDEResult

Returns:

SDEResult

Examples

>>> import numpy as np
>>> result = geometric_brownian_motion(0.1, 0.4, n_steps=50, n_paths=40000, seed=0)
>>> bool(abs(result.paths[:, -1].mean() - np.exp(0.1)) < 0.01)  # E X_t = X_0 e^{mu t}
True
mathematicskit.probability.gibbs_sampler(conditional_samplers, x0, n_samples=10000, burn_in=0, seed=0)[source]#

Sample a joint distribution by redrawing each coordinate from its full conditional.

One sweep (systematic scan) updates \(x_1, \dots, x_d\) in turn, drawing \(x_i \sim \pi(x_i \mid x_{-i})\) with the other coordinates at their current values. Each update leaves \(\pi\) invariant, so the sweep does too (Geman and Geman, 1984).

Parameters:

  • conditional_samplers (Sequence[Callable[[ndarray, Generator], float]]) – conditional_samplers[i](x, rng) -> float draws coordinate i from its conditional distribution given the rest of x.

  • x0 (array_like, shape (d,)) – Starting state, with d == len(conditional_samplers).

  • n_samples (int) – Number of sweeps kept after burn-in.

  • burn_in (int) – Number of initial sweeps discarded.

  • seed (int)

Return type:

MCMCResult

Returns:

MCMCResult – One row per sweep; acceptance_rate is 1.

Examples

>>> import numpy as np
>>> rho = 0.8  # bivariate normal: x | y ~ N(rho y, 1 - rho^2), and symmetrically
>>> s = np.sqrt(1 - rho**2)
>>> conditionals = [lambda x, rng: rho * x[1] + s * rng.standard_normal(), lambda x, rng: rho * x[0] + s * rng.standard_normal()]
>>> result = gibbs_sampler(conditionals, [0.0, 0.0], n_samples=40000, seed=0)
>>> bool(abs(np.corrcoef(result.samples.T)[0, 1] - rho) < 0.02)
True
mathematicskit.probability.importance_sampling_integrate(f, proposal_sampler, proposal_pdf, n=10000, seed=0)[source]#

Importance-sampling estimate of \(\int f(x)\,dx = E_q[f(X)/q(X)]\) for \(X \sim q\).

Draws from an arbitrary proposal density \(q\) (via proposal_sampler) rather than uniformly, reweighting each sample by \(f(x)/q(x)\); when \(q\) is chosen proportional to \(|f|\), this can drastically reduce variance versus plain Monte Carlo (monte_carlo_integrate()) for integrands concentrated in a small region. See Robert & Casella, Monte Carlo Statistical Methods, 2nd ed., Ch. 3.3.

Parameters:
Return type:

MonteCarloResult

Returns:

MonteCarloResult

Examples

>>> import numpy as np
>>> from mathematicskit.probability.systems.continuous import Exponential
>>> # Integrate exp(-x^2/2) over [0, inf) using an Exponential(1) proposal
>>> # (concentrated where the integrand is largest, near 0).
>>> proposal = Exponential(rate=1.0)
>>> f = lambda x: np.exp(-x**2 / 2.0)
>>> result = importance_sampling_integrate(f, lambda rng, n: proposal.sample(size=n, seed=rng), proposal.pdf, n=200000, seed=0)
>>> bool(abs(result.estimate - np.sqrt(np.pi / 2.0)) < 0.02)
True
mathematicskit.probability.integrated_autocorrelation_time(chain, window_factor=5.0)[source]#

Integrated autocorrelation time \(\tau = 1 + 2\sum_{k \geq 1} \rho_k\) of a Markov chain.

The variance of a mean over \(n\) correlated draws is \(\tau\) times that over \(n\) independent ones. The sum is truncated at the smallest window \(M \geq c\,\hat\tau(M)\) (Sokal’s automatic windowing, \(c\) = window_factor), which balances the bias of stopping early against the noise of the far-lag estimates. The autocorrelation is computed with numpy.fft.rfft(). See A. D. Sokal, “Monte Carlo Methods in Statistical Mechanics: Foundations and New Algorithms,” in Functional Integration (Springer, 1997), Sec. 3.

Parameters:

  • chain (ndarray) – A scalar function of the chain’s states, in order.

  • window_factor (float)

Return type:

float

Returns:

float

Examples

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> bool(abs(integrated_autocorrelation_time(rng.standard_normal(100000)) - 1.0) < 0.05)
True
>>> # AR(1) with coefficient phi has tau = (1 + phi) / (1 - phi) = 9 for phi = 0.8.
>>> from scipy.signal import lfilter
>>> x = lfilter([1.0], [1.0, -0.8], rng.standard_normal(200000))
>>> bool(abs(integrated_autocorrelation_time(x) - 9.0) < 0.5)
True
mathematicskit.probability.law_of_large_numbers_trace(dist, n_values, seed=0)[source]#

Sample means at increasing sample sizes, illustrating the Law of Large Numbers.

Draws one long sequence of max(n_values) i.i.d. samples from dist and returns the running sample mean truncated at each requested n – by the (strong) LLN, this sequence converges almost surely to \(E[X]\) as \(n \to \infty\). See DeGroot & Schervish, Probability and Statistics, 4th ed., Theorem 6.2.4.

Parameters:
Return type:

ndarray

Returns:

ndarray, shape (len(n_values),) – The sample mean of the first n draws, for each n in n_values.

Examples

>>> from mathematicskit.probability.systems.continuous import Exponential
>>> means = law_of_large_numbers_trace(Exponential(rate=2.0), n_values=[100, 10000, 1000000], seed=0)
>>> bool(abs(means[-1] - 0.5) < abs(means[0] - 0.5))
True
mathematicskit.probability.metropolis_hastings(log_target, x0, n_samples=10000, proposal_scale=1.0, burn_in=0, seed=0, proposal=None, log_proposal_density=None)[source]#

Sample from \(\pi \propto e^{\text{log\_target}}\) with the Metropolis-Hastings algorithm.

By default the proposal is a Gaussian random walk \(y = x + \sigma Z\), which is symmetric, so the acceptance ratio reduces to Metropolis et al.’s (1953) \(\pi(y)/\pi(x)\). Passing proposal and log_proposal_density gives Hastings’s (1970) general, possibly asymmetric, version.

Parameters:

  • log_target (Callable[[ndarray], float]) – log_target(x) -> float, the log of the target density up to an additive constant; -inf outside its support.

  • x0 (array_like) – Starting state, with log_target(x0) finite.

  • n_samples (int) – Number of states kept after burn-in.

  • proposal_scale (float) – Standard deviation \(\sigma\) of the default random-walk proposal.

  • burn_in (int) – Number of initial steps discarded.

  • seed (int)

  • proposal (Callable[[ndarray, Generator], ndarray] | None) – proposal(x, rng) -> y, a draw from \(q(\cdot \mid x)\).

  • log_proposal_density (Callable[[ndarray, ndarray], float] | None) – log_proposal_density(y, x) -> float, \(\log q(y \mid x)\) up to a constant. Required with an asymmetric proposal; if omitted, the proposal is taken to be symmetric.

Return type:

MCMCResult

Returns:

MCMCResult

Examples

>>> import numpy as np
>>> result = metropolis_hastings(lambda x: -0.5 * float(x @ x), [3.0], n_samples=40000, proposal_scale=2.4, seed=0)
>>> bool(abs(result.samples.mean()) < 0.05), bool(abs(result.samples.var() - 1.0) < 0.05)
(True, True)
mathematicskit.probability.monte_carlo_integrate(f, a, b, n=10000, seed=0)[source]#

Plain Monte Carlo estimate of \(\int_a^b f(x)\,dx\).

\(\int_a^b f\,dx = (b-a)\,E_{X\sim U(a,b)}[f(X)]\), estimated by the sample mean of \(f\) at n points drawn uniformly from [a, b]; the standard error follows from the sample variance via the CLT. See Robert & Casella, Monte Carlo Statistical Methods, 2nd ed., Ch. 3.2.

Parameters:
Return type:

MonteCarloResult

Returns:

MonteCarloResult

Examples

>>> result = monte_carlo_integrate(lambda x: x**2, 0.0, 1.0, n=200000, seed=0)
>>> abs(result.estimate - 1.0 / 3.0) < 0.01
True
mathematicskit.probability.ornstein_uhlenbeck(theta, mu, sigma, x0=0.0, t_max=1.0, n_steps=1000, n_paths=1, seed=0)[source]#

Sample the Ornstein-Uhlenbeck process \(dX = \theta(\mu - X)\,dt + \sigma\,dW\) by Euler-Maruyama.

The process reverts to \(\mu\) at rate \(\theta\). It is Gaussian, with

\[E X_t = \mu + (X_0 - \mu)e^{-\theta t}, \qquad \operatorname{Var} X_t = \frac{\sigma^2}{2\theta}\left(1 - e^{-2\theta t}\right),\]

and stationary distribution \(\mathcal N(\mu, \sigma^2/2\theta)\) (Uhlenbeck and Ornstein, 1930).

Parameters:

  • theta (float) – Mean-reversion rate, theta > 0.

  • mu (float) – Long-run mean.

  • sigma (float) – Noise amplitude.

  • x0 (float)

  • t_max (float)

  • n_steps (int)

  • n_paths (int)

  • seed (int)

Return type:

SDEResult

Returns:

SDEResult

Examples

>>> result = ornstein_uhlenbeck(2.0, 1.0, 0.5, x0=3.0, t_max=5.0, n_steps=500, n_paths=20000)
>>> bool(abs(result.paths[:, -1].mean() - 1.0) < 0.01), bool(abs(result.paths[:, -1].var() - 0.0625) < 0.005)
(True, True)
mathematicskit.probability.random_walk_return_fraction(dim, n_steps, n_walks=2000, seed=0)[source]#

Fraction of simple random walks on \(\mathbb{Z}^d\) that revisit the origin within n_steps steps.

Estimates Pólya’s return probability truncated at n_steps. The walks are advanced one step at a time, so memory stays \(O(n_\text{walks}\,d)\) however long they run.

Parameters:
Return type:

float

Returns:

float

Examples

>>> frac = random_walk_return_fraction(dim=1, n_steps=100, n_walks=20000, seed=0)
>>> bool(abs(frac - return_probability_1d(100)) < 0.01)
True
mathematicskit.probability.return_probability_1d(n_steps)[source]#

Exact probability that a 1D simple random walk revisits 0 within n_steps steps.

A walk avoids 0 for its first \(2m\) steps with probability \(\binom{2m}{m}4^{-m} \sim 1/\sqrt{\pi m}\), so the return probability is

\[P(\text{return by } 2m) = 1 - \binom{2m}{m}4^{-m} \to 1,\]

Pólya’s recurrence in one dimension (Feller, An Introduction to Probability Theory and Its Applications, vol. 1, 3rd ed., Sec. III.3).

Parameters:

n_steps (int) – Number of steps; an odd count gives the same answer as n_steps - 1.

Return type:

float

Returns:

float

Examples

>>> return_probability_1d(2)
0.5
>>> round(return_probability_1d(4), 4)
0.625
mathematicskit.probability.rule_of_succession(successes, trials)[source]#

Laplace’s rule of succession, \(P(\text{next success}) = \dfrac{k+1}{n+2}\).

The posterior predictive probability of a success on trial \(n + 1\) under Bayes’s uniform prior: the mean of beta_binomial_posterior() (k, n). Laplace (1774) famously applied it to the probability that the sun will rise tomorrow.

Parameters:

  • successes (int)

  • trials (int)

Return type:

float

Returns:

float

Examples

>>> rule_of_succession(0, 0)
0.5
>>> round(rule_of_succession(9, 10), 4)
0.8333
mathematicskit.probability.simple_random_walk(n_walks=1, n_steps=1000, dim=1, seed=0)[source]#

Positions of simple random walks on \(\mathbb{Z}^d\) started at the origin.

Parameters:
Return type:

ndarray

Returns:

ndarray of int, shape (n_walks, n_steps + 1, dim)

Examples

>>> walks = simple_random_walk(n_walks=3, n_steps=10, dim=2, seed=0)
>>> walks.shape
(3, 11, 2)
>>> bool(np.all(np.abs(np.diff(walks, axis=1)).sum(axis=2) == 1))
True
mathematicskit.probability.st_petersburg_certainty_equivalent(wealth=0.0, n_terms=200)[source]#

The sure amount a log-utility player values the St. Petersburg game at.

With initial wealth \(w\), Bernoulli’s player is indifferent between playing and receiving the sure amount \(c\) solving

\[\ln(w + c) = \sum_{k=1}^{\infty} 2^{-k} \ln(w + 2^k).\]

For \(w = 0\) the right side is \(2\ln 2\), so \(c = 4\): a game with infinite expected value is worth just four ducats. The series is truncated after n_terms terms, far past double precision for any reasonable wealth.

Parameters:

  • wealth (float) – The player’s initial wealth \(w \geq 0\).

  • n_terms (int) – Number of series terms.

Return type:

float

Returns:

float – The certainty equivalent \(c\).

Examples

>>> round(st_petersburg_certainty_equivalent(0.0), 10)
4.0
mathematicskit.probability.st_petersburg_payoffs(n_games, seed=0)[source]#

Simulated payoffs \(2^K\) of the St. Petersburg game, with \(K \sim \mathrm{Geometric}(1/2)\).

Parameters:

  • n_games (int) – Number of independent games played.

  • seed (int)

Return type:

ndarray

Returns:

ndarray, shape (n_games,) – The payoff of each game, as floats.

Examples

>>> payoffs = st_petersburg_payoffs(5, seed=0)
>>> bool(payoffs.min() >= 2.0)
True
On this page
Show Source