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, and continuous-time Markov chains via scipy.linalg.expm().

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.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.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.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.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.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.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.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.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