Bayesian Statistics

bayesian-statistics-ai999x

Bayes' theorem, conjugate priors and posteriors, credible intervals, Bayes factors, Monte Carlo estimation, uncertainty propagation, sensitivity analysis and decision analysis.

Mathematics & Statistics Question Answering Sentence Similarity Feature Extraction

Hugging Face ↗GitHub ↗See it on the capability orb

78capabilities

12topics

78formulas typeset

8worked examples below

1.0.0version

About this model

What Bayesian Statistics does

Bayesian Statistics (bayesian-statistics-ai999x) is an offline Python package for bayesian statistics. It covers Bayes' theorem, conjugate priors and posteriors, credible intervals, Bayes factors, Monte Carlo estimation, uncertainty propagation, sensitivity analysis and decision analysis.

It has 78 capabilities — 78 computational procedures — in 12 topics, the largest being Conjugate priors, Uncertainty propagation, Decision analysis, Credible intervals, Bayes factors and model comparison and Monte Carlo estimation. Ask a question in words with ask(), or call a capability by name with typed inputs with solve(). Every answer carries its value and unit, the steps that produced it and the formula in LaTeX.

When a question is not determined — a value is missing, an input is out of range, or the request is outside the package — it says so and names what it needs, instead of guessing. It runs on an ordinary CPU without a network connection (Python >=3.10), and the same question always gets the same answer.

Install

pip install bayesian-statistics-ai999x
hf download AI999X/bayesian-statistics-ai999x --local-dir bayesian-statistics-ai999x

What it covers

Bayes factors and model comparison6Bayes factor for a normal mean (known variance) · Bayes factor for a point null on a proportion (beta prior under H1) · Bayes factor for two simple binomial hypotheses …Bayes' theorem and conditional probability5Bayes' theorem for an event and its complement · Law of total probability · Positive and negative predictive values of a test …Conjugate priors17Bayesian A/B test: probability that A beats B · Beta-binomial conjugate posterior · Dirichlet-multinomial conjugate posterior …Credible intervals7Credible interval for a Poisson rate · Credible interval for a proportion (beta-binomial) · Equal-tailed credible interval of a beta distribution …Decision analysis9Bayes action: minimum posterior expected loss · Best action by expected monetary value · Break-even (indifference) probability of two actions …Information and quadrature for Bayesian computation3Gauss-Hermite rule for expectations under a normal distribution · Kullback-Leibler divergence between two beta distributions · Kullback-Leibler divergence between two normal distributions Monte Carlo estimation6Control-variate Monte Carlo estimate · Importance sampling estimate and effective sample size · Monte Carlo estimate of an integral from uniform samples …Odds and likelihood ratios4Likelihood ratios of a diagnostic test · Odds and log-odds of a probability · Post-test probability from a likelihood ratio …Posterior approximation1Laplace (normal) approximation of a beta posterior Prior elicitation1Beta prior from a mean and standard deviation Sensitivity analysis5Morris elementary-effects screening · Sobol indices of a linear (additive) model · Sobol indices of a model with an interaction term …Uncertainty propagation14Covariance propagation through a linear map y = A x · Effective degrees of freedom and expanded uncertainty (Welch-Satterthwaite) · Exact mean and standard deviation of a product of independent normals …

Try it

Worked examples

Each question below was put to bayesian-statistics-ai999x and the answer shown is the one it returned, with the steps it took.

Bayes' theorem and conditional probabilitym.ask(…)

A test has sensitivity 0.9 and specificity 0.95; the prevalence is 0.01. What is the positive predictive value?

ppv = 0.1538461538

npv = 0.9989378651

p_positive = 0.0585

false_discovery_rate = 0.8461538462

false_omission_rate = 0.001062134891

Show the steps

Positive and negative predictive values of a test
\mathrm{PPV} = \frac{\mathrm{Se}\,\pi}{\mathrm{Se}\,\pi + (1-\mathrm{Sp})(1-\pi)},\quad \mathrm{NPV} = \frac{\mathrm{Sp}\,(1-\pi)}{\mathrm{Sp}\,(1-\pi) + (1-\mathrm{Se})\,\pi}
with prevalence = 0.01, sensitivity = 0.9, specificity = 0.95
→ ppv = 0.1538461538, npv = 0.9989378651, p_positive = 0.0585, false_discovery_rate = 0.8461538462, false_omission_rate = 0.001062134891

Capability: Positive and negative predictive values of a test

Conjugate priorsm.ask(…)

Beta-binomial posterior: prior Beta(2, 2), 7 successes in 10 trials. Find the posterior mean.

post_alpha = 9

post_beta = 5

posterior_mean = 0.6428571429

posterior_variance = 0.01530612245

posterior_sd = 0.1237179148

Show the steps

Beta-binomial conjugate posterior
p\mid k \sim \mathrm{Beta}(\alpha + k,\ \beta + n - k),\quad E[p\mid k] = \frac{\alpha + k}{\alpha + \beta + n}
with alpha = 2, beta = 2, successes = 7, trials = 10
→ post_alpha = 9, post_beta = 5, posterior_mean = 0.6428571429, posterior_variance = 0.01530612245, posterior_sd = 0.1237179148

Capability: Beta-binomial conjugate posterior

Credible intervalsm.ask(…)

95% highest-density interval of a Beta(3, 9) posterior.

lower = 0.04055516641

upper = 0.4837236593

width = 0.4431684929

Show the steps

Highest-density interval of a beta posterior
F(u) - F(l) = c,\quad f(l) = f(u)
with alpha = 3, beta = 9, level = 0.95
→ lower = 0.04055516641, upper = 0.4837236593, width = 0.4431684929

Capability: Highest-density interval of a beta posterior

Decision analysism.ask(…)

Expected value of perfect information: payoffs = [[100, -20], [40, 30]], probabilities = [0.4, 0.6].

evpi = 24

ev_with_perfect_information = 58

best_expected_value = 34

Show the steps

Expected value of perfect information (EVPI)
EVPI = \sum_s P(s)\max_a u(a,s) - \max_a\sum_s P(s)\,u(a,s)
with payoffs = [100 -20; 40 30], probabilities = [0.4, 0.6]
→ evpi = 24, ev_with_perfect_information = 58, best_expected_value = 34

Capability: Expected value of perfect information (EVPI)

Uncertainty propagationm.solve(…)

m.solve("uncertainty_product_quotient", a=12.0, ua=0.3, b=4.0, ub=0.2)

product = 48

u_product = 2.683281573

quotient = 3

u_quotient = 0.1677050983

Show the steps

Uncertainty of a product and a quotient
\frac{u_{ab}}{|ab|} = \frac{u_{a/b}}{|a/b|} = \sqrt{\left(\frac{u_a}{a}\right)^2 + \left(\frac{u_b}{b}\right)^2}
with a = 12, ua = 0.3, b = 4, ub = 0.2
→ product = 48, u_product = 2.683281573, quotient = 3, u_quotient = 0.1677050983

Capability: Uncertainty of a product and a quotient

Bayes factors and model comparisonm.ask(…)

Bayes factor for p1 = 0.7 against p0 = 0.5 after 14 successes in 20 trials.

bf10 = 5.184417905

log10_bf10 = 0.7147000018

posterior_h1 = 0.8383032946

Show the steps

Bayes factor for two simple binomial hypotheses
BF_{10} = \left(\frac{p_1}{p_0}\right)^{k}\left(\frac{1-p_1}{1-p_0}\right)^{n-k}
with p0 = 0.5, p1 = 0.7, successes = 14, trials = 20
→ bf10 = 5.184417905, log10_bf10 = 0.7147000018, posterior_h1 = 0.8383032946

Capability: Bayes factor for two simple binomial hypotheses

Information and quadrature for Bayesian computationm.ask(…)

Gauss-Hermite rule with n = 3 nodes for a normal expectation.

nodes = [-1.732050808, 0, 1.732050808]

weights = [0.1666666667, 0.6666666667, 0.1666666667]

Show the steps

Gauss-Hermite rule for expectations under a normal distribution
E[f(Z)] \approx \sum_{i=1}^{n} w_i f(x_i),\quad He_n(x_i) = 0,\quad w_i = \frac{n!}{n^2\,He_{n-1}(x_i)^2}
with n = 3
→ nodes = [-1.732050808, 0, 1.732050808], weights = [0.1666666667, 0.6666666667, 0.1666666667]

Capability: Gauss-Hermite rule for expectations under a normal distribution

Monte Carlo estimationm.ask(…)

Control variate: outputs = [1, 2, 3], controls = [1, 2, 3], known mean = 2.5.

estimate = 2.5

coefficient = 1

plain_estimate = 2

variance_reduction = 0

Show the steps

Control-variate Monte Carlo estimate
\hat\mu_{cv} = \bar y - \hat b\,(\bar c - \mu_c),\quad \hat b = \frac{\sum (y_i-\bar y)(c_i-\bar c)}{\sum (c_i-\bar c)^2}
with y = [1, 2, 3], c = [1, 2, 3], control_mean = 2.5
→ estimate = 2.5, coefficient = 1, plain_estimate = 2, variance_reduction = 0

Capability: Control-variate Monte Carlo estimate

Use it

From Python or the command line

Every snippet was run as shown; the output under it is what it printed.

import bayesian_statistics_ai999x as m

a = m.ask("A test has sensitivity 0.9 and specificity 0.95; the prevalence is 0.01. What is the positive predictive value?")
print(a.outcome)
print(a.values)
for step in a.steps:
    print(step["name"])

Output

SOLVED
{'ppv': 0.153846153846, 'npv': 0.998937865109, 'p_positive': 0.0585, 'false_discovery_rate': 0.846153846154, 'false_omission_rate': 0.00106213489113}
Positive and negative predictive values of a test

import bayesian_statistics_ai999x as m

a = m.solve("uncertainty_product_quotient", a=12.0, ua=0.3, b=4.0, ub=0.2)
print(a.outcome)
print(a.values)

Output

SOLVED
{'product': 48.0, 'u_product': 2.683281573, 'quotient': 3.0, 'u_quotient': 0.167705098312}

import bayesian_statistics_ai999x as m

# an input is missing: the answer says what it needs instead of guessing
a = m.solve("uncertainty_product_quotient")
print(a.outcome)
print(a.message)

Output

NEEDS_INPUT
uncertainty_product_quotient needs a, ua, b, ub (value of a; standard uncertainty of a (>= 0); value of b (non-zero); standard uncertainty of b (>= 0)). Example: solve('uncertainty_product_quotient', a=12, ua=0.3, b=4, ub=0.2).

bayesian-statistics-ai999x "A test has sensitivity 0.9 and specificity 0.95; the prevalence is 0.01. What is the positive predictive value?"

Output

ppv = 0.153846, npv = 0.998938, p_positive = 0.0585, false_discovery_rate = 0.846154, false_omission_rate = 0.00106213
1. Positive and negative predictive values of a test: Bayes' theorem with the total probability of a positive result  with prevalence = 0.01, sensitivity = 0.9, specificity = 0.95  ->  ppv = 0.153846, npv = 0.998938, p_positive = 0.0585, false_discovery_rate = 0.846154, false_omission_rate = 0.00106213
Assumptions: probabilities between 0 and 1; both results possible
…

Figures

Computed, not drawn by hand

Every value in these figures was computed by the package, each point a run of solve(); in a notebook, plot() and plot_relation() draw a capability's curves the same way.

Beta prior and posterior after 7 successes in 10 trials

p\mid k \sim \mathrm{Beta}(\alpha + k,\ \beta + n - k),\quad E[p\mid k] = \frac{\alpha + k}{\alpha + \beta + n}

A Beta(2, 2) prior updated with 7 successes in 10 trials: the posterior Beta(9, 5) with mean 0.642857 and standard deviation 0.123718; shaded: the equal-tailed 95 % credible interval [0.385738, 0.861421]. Both densities are drawn from the package's cumulative probabilities at 199 values of θ (each interval's probability divided by its width).

Posterior mean for every count of successes and failures

p\mid k \sim \mathrm{Beta}(\alpha + k,\ \beta + n - k),\quad E[p\mid k] = \frac{\alpha + k}{\alpha + \beta + n}

The posterior mean of a success probability under a Beta(2, 2) prior after s successes (rows, 0 to 5) and f failures (columns, 0 to 5), each cell from the package (rounded to 3 decimals): from 0.222 after 5 failures to 0.778 after 5 successes, 0.5 with no data.

Probability of the condition after repeated tests

\text{odds}_{post} = \frac{\pi}{1-\pi}\,(LR^{+})^{k}\,(LR^{-})^{m}

Repeated diagnostic tests with prevalence 0.01, sensitivity 0.9 and specificity 0.95: the package's posterior probability of the condition after k = 1 … 5 positive results (0.153846, 0.765957, 0.983308, 0.999058, 0.999948) and after 1 … 5 negative results (down to 1.3054e-07).

Prior and posterior over three hypotheses

P(H_i\mid D) = \frac{P(D\mid H_i)\,P(H_i)}{\sum_j P(D\mid H_j)\,P(H_j)}

Bayes' rule over 3 hypotheses with priors 0.5, 0.3, 0.2 and likelihoods 0.1, 0.4, 0.7: the package's posterior probabilities 0.16129, 0.387097, 0.451613 (evidence 0.31); the most probable hypothesis is 3 with 0.451613.

Category probabilities after multinomial counts

\theta\mid n \sim \mathrm{Dir}(\alpha_1 + n_1, \ldots, \alpha_k + n_k),\quad E[\theta_j] = \frac{\alpha_j + n_j}{\alpha_0 + n}

A Dirichlet(1, 1, 1) prior over 3 categories updated with the counts 10, 5, 3: the package's posterior Dirichlet(11, 6, 4), posterior means 0.52381, 0.285714, 0.190476; whiskers ±2 posterior standard deviations (0.106479, 0.0963143, 0.083719).

Normal posterior shrinking towards the data as n grows

\mu_n = B\mu_0 + (1-B)\bar x,\ B = \frac{\sigma^2/n}{\sigma^2/n + \tau_0^2},\quad \mu_n \pm z_{(1+c)/2}\,\tau_n

A normal prior (mean 0, sd 1) and data with σ = 2 and sample mean 1: the package's posterior mean and 95 % credible interval for 17 sample sizes from 1 to 100, moving from the prior towards the data; marked: n = 10, posterior mean = 0.714286, interval [-0.333359, 1.76193].

Formulas

The mathematics inside

12 of its 78 formulas; every one is in the capability browser below.

Positive and negative predictive values of a test

\mathrm{PPV} = \frac{\mathrm{Se}\,\pi}{\mathrm{Se}\,\pi + (1-\mathrm{Sp})(1-\pi)},\quad \mathrm{NPV} = \frac{\mathrm{Sp}\,(1-\pi)}{\mathrm{Sp}\,(1-\pi) + (1-\mathrm{Se})\,\pi}

Beta-binomial conjugate posterior

p\mid k \sim \mathrm{Beta}(\alpha + k,\ \beta + n - k),\quad E[p\mid k] = \frac{\alpha + k}{\alpha + \beta + n}

Highest-density interval of a beta posterior

F(u) - F(l) = c,\quad f(l) = f(u)

Expected value of perfect information (EVPI)

EVPI = \sum_s P(s)\max_a u(a,s) - \max_a\sum_s P(s)\,u(a,s)

Uncertainty of a product and a quotient

\frac{u_{ab}}{|ab|} = \frac{u_{a/b}}{|a/b|} = \sqrt{\left(\frac{u_a}{a}\right)^2 + \left(\frac{u_b}{b}\right)^2}

Bayes factor for two simple binomial hypotheses

BF_{10} = \left(\frac{p_1}{p_0}\right)^{k}\left(\frac{1-p_1}{1-p_0}\right)^{n-k}

Gauss-Hermite rule for expectations under a normal distribution

E[f(Z)] \approx \sum_{i=1}^{n} w_i f(x_i),\quad He_n(x_i) = 0,\quad w_i = \frac{n!}{n^2\,He_{n-1}(x_i)^2}

Control-variate Monte Carlo estimate

\hat\mu_{cv} = \bar y - \hat b\,(\bar c - \mu_c),\quad \hat b = \frac{\sum (y_i-\bar y)(c_i-\bar c)}{\sum (c_i-\bar c)^2}

Bayes factor for a normal mean (known variance)

BF_{01} = \frac{p(\mu_0\mid \bar x, H_1)}{p(\mu_0\mid H_1)}

Bayes' theorem for an event and its complement

P(A\mid B) = \frac{P(B\mid A)\,P(A)}{P(B\mid A)\,P(A) + P(B\mid \neg A)\,(1 - P(A))}

Bayesian A/B test: probability that A beats B

P(p_B > p_A) = \sum_{i=0}^{\alpha_B - 1}\frac{B(\alpha_A + i,\ \beta_A + \beta_B)}{(\beta_B + i)\,B(1 + i,\ \beta_B)\,B(\alpha_A, \beta_A)}

Credible interval for a Poisson rate

\left[\frac{G^{-1}_{\alpha'}\!\left(\frac{1-c}{2}\right)}{\beta'},\ \frac{G^{-1}_{\alpha'}\!\left(\frac{1+c}{2}\right)}{\beta'}\right]

Capabilities

Every capability, with an example

Search by name, formula, input or topic. Each capability opens to show what it takes, what it returns, when it applies, and a question it answered.

78 capabilities

Bayes factors and model comparison 6

Bayes factor for a normal mean (known variance)bayes_factor_normal_mean

BF10 for H1: mu ~ N(mu0, tau^2) against H0: mu = mu0, from the sample mean of n observations with known sigma; BF01 and the posterior probability of H1 at equal prior odds.

BF_{01} = \frac{p(\mu_0\mid \bar x, H_1)}{p(\mu_0\mid H_1)}

Inputs

sample_mean sample mean
n sample size (>= 1) integer
sigma known standard deviation (> 0)
mu0 null value mu0
prior_sd prior standard deviation tau of mu under H1 (> 0)

Returns

bf10 Bayes factor H1 : H0
bf01 1/BF10
log10_bf10 log10 BF10
posterior_h1 P(H1 | data) with prior odds 1

Applies when

sigma, tau > 0
n >= 1

Example

m.ask("Bayes factor for a normal mean: sample mean = 0, n = 1, sigma = 1, mu0 = 0, prior sd = 1.")

→ bf10 = 0.707107, bf01 = 1.41421, log10_bf10 = -0.150515, posterior_h1 = 0.414214

Bayes factor for a point null on a proportion (beta prior under H1)bayes_factor_point_null_binomial

BF10 for H1: p ~ Beta(alpha, beta) against H0: p = p0 after k successes in n trials (default uniform prior), with BF01 and the posterior probability of H1 at equal prior odds.

BF_{10} = \frac{B(\alpha + k,\ \beta + n - k)/B(\alpha,\beta)}{p_0^{k}(1-p_0)^{n-k}}

Inputs

successes successes k integer
trials trials n (>= 1) integer
p0 null value p0
alpha Beta alpha under H1 (default 1) optional
beta Beta beta under H1 (default 1) optional

Returns

bf10 Bayes factor H1 : H0
bf01 1/BF10
log10_bf10 log10 BF10
posterior_h1 P(H1 | data) with prior odds 1

Applies when

0 < p0 < 1
alpha, beta > 0
0 <= k <= n

Example

m.ask("Bayes factor against the point null p0 = 0.5 with a uniform prior: 2 successes in 2 trials.")

→ bf10 = 1.33333, bf01 = 0.75, log10_bf10 = 0.124939, posterior_h1 = 0.571429

Bayes factor for two simple binomial hypothesesbayes_factor_binomial_simple

Bayes factor BF10 for H1: p = p1 against H0: p = p0 after k successes in n trials, its log10 and the posterior probability of H1 with equal prior odds.

BF_{10} = \left(\frac{p_1}{p_0}\right)^{k}\left(\frac{1-p_1}{1-p_0}\right)^{n-k}

Inputs

successes successes k integer
trials trials n (>= 1) integer
p0 success probability under H0
p1 success probability under H1

Returns

bf10 Bayes factor H1 : H0
log10_bf10 log10 of BF10
posterior_h1 P(H1 | data) with prior odds 1

Applies when

0 < p0, p1 < 1
0 <= k <= n

Example

m.ask("Bayes factor for p1 = 0.7 against p0 = 0.5 after 14 successes in 20 trials.")

→ bf10 = 5.18442, log10_bf10 = 0.7147, posterior_h1 = 0.838303

Marginal likelihood (evidence) of binomial data under a beta priorbeta_binomial_evidence

Prior predictive probability of k successes in n trials when p ~ Beta(alpha, beta), and its logarithm.

p(k) = \binom{n}{k}\frac{B(\alpha + k,\ \beta + n - k)}{B(\alpha,\beta)}

Inputs

alpha prior Beta alpha (> 0)
beta prior Beta beta (> 0)
successes number of successes k integer
trials number of trials n integer

Returns

marginal_likelihood P(k | n, alpha, beta)
log_marginal_likelihood its natural log

Applies when

alpha, beta > 0
0 <= k <= n

Example

m.ask("Marginal likelihood of 7 successes in 10 trials under a Beta(1, 1) prior.")

→ marginal_likelihood = 0.0909091, log_marginal_likelihood = -2.3979

Posterior model probabilities from log evidencesposterior_model_probabilities

Posterior probabilities of several models from their log marginal likelihoods and prior model probabilities (default equal), computed stably in log space.

P(M_i\mid D) = \frac{P(D\mid M_i)P(M_i)}{\sum_j P(D\mid M_j)P(M_j)}

Inputs

log_evidences natural-log marginal likelihoods of the models vector
priors prior model probabilities (default equal) vector optional

Returns

posterior posterior model probabilities
best_model number of the most probable model
best_probability its posterior probability

Applies when

prior probabilities sum to 1

Example

m.ask("Posterior model probabilities from log marginal likelihoods = [-10, -10, -10].")

→ posterior = [0.333333333333, 0.333333333333, 0.333333333333], best_model = 1, best_probability = 0.333333

Posterior probability of a hypothesis from a Bayes factorposterior_from_bayes_factor

Posterior probability of H1 from the Bayes factor BF10 and the prior probability of H1 (default 1/2).

\frac{P(H_1\mid D)}{P(H_0\mid D)} = BF_{10}\,\frac{P(H_1)}{P(H_0)}

Inputs

bayes_factor Bayes factor BF10 in favour of H1 (> 0)
prior prior probability of H1 (default 0.5) optional

Returns

posterior posterior probability of H1
posterior_odds posterior odds H1 : H0
posterior_h0 posterior probability of H0

Applies when

BF10 > 0
0 < prior < 1

Example

m.ask("A Bayes factor of 6 with prior probability 0.5: what is the posterior probability of H1?")

→ posterior = 0.857143, posterior_odds = 6, posterior_h0 = 0.142857

Bayes' theorem and conditional probability 5

Bayes' theorem for an event and its complementbayes_theorem

The posterior probability P(A|B) from the prior P(A) and the likelihoods P(B|A) and P(B|not A).

P(A\mid B) = \frac{P(B\mid A)\,P(A)}{P(B\mid A)\,P(A) + P(B\mid \neg A)\,(1 - P(A))}

Inputs

p_a prior probability P(A)
p_b_given_a likelihood P(B|A)
p_b_given_not_a likelihood P(B|not A)

Returns

p_a_given_b posterior probability P(A|B)
p_b total probability P(B)
p_not_a_given_b P(not A|B) = 1 - P(A|B)

Applies when

probabilities between 0 and 1
P(B) > 0

Example

m.ask("Bayes' theorem: P(A) = 0.3, P(B|A) = 0.8, P(B|not A) = 0.1. Find P(A|B).")

→ p_a_given_b = 0.774194, p_b = 0.31, p_not_a_given_b = 0.225806

Law of total probabilitytotal_probability

P(B) = sum of P(B|A_i) P(A_i) over a partition A_1, ..., A_k (and the joint probabilities).

P(B) = \sum_i P(B\mid A_i)\,P(A_i)

Inputs

priors probabilities P(A_i) of the partition (sum to 1) vector
conditionals conditional probabilities P(B|A_i), one per part vector

Returns

p_b total probability P(B)
joint joint probabilities P(B and A_i)

Applies when

the P(A_i) sum to 1
same number of conditionals as parts

Example

m.ask("Law of total probability: priors = [0.5, 0.3, 0.2], conditionals = [0.01, 0.02, 0.05]. Find P(B).")

→ p_b = 0.021, joint = [0.005, 0.006, 0.01]

Positive and negative predictive values of a testdiagnostic_test_predictive_values

The probability of the condition after a positive result (PPV) and of no condition after a negative one (NPV), from prevalence, sensitivity and specificity.

\mathrm{PPV} = \frac{\mathrm{Se}\,\pi}{\mathrm{Se}\,\pi + (1-\mathrm{Sp})(1-\pi)},\quad \mathrm{NPV} = \frac{\mathrm{Sp}\,(1-\pi)}{\mathrm{Sp}\,(1-\pi) + (1-\mathrm{Se})\,\pi}

Inputs

prevalence prevalence (pre-test probability of the condition)
sensitivity sensitivity P(positive | condition)
specificity specificity P(negative | no condition)

Returns

ppv positive predictive value P(condition | positive)
npv negative predictive value P(no condition | negative)
p_positive probability of a positive result
false_discovery_rate 1 - PPV
false_omission_rate 1 - NPV

Applies when

probabilities between 0 and 1
both results possible

Example

m.ask("A test has sensitivity 0.9 and specificity 0.95; the prevalence is 0.01. What is the positive predictive value?")

→ ppv = 0.153846, npv = 0.998938, p_positive = 0.0585, false_discovery_rate = 0.846154, false_omission_rate = 0.00106213

Posterior after repeated independent testssequential_test_posterior

Probability of the condition after several conditionally independent test results (k positive and m negative), by multiplying the pre-test odds by LR+ for each positive and LR- for each negative.

\text{odds}_{post} = \frac{\pi}{1-\pi}\,(LR^{+})^{k}\,(LR^{-})^{m}

Inputs

prevalence prevalence (pre-test probability of the condition)
sensitivity sensitivity P(positive | condition)
specificity specificity P(negative | no condition)
positives number of positive results integer
negatives number of negative results (default 0) integer optional

Returns

posterior probability of the condition after all results
posterior_odds its odds

Applies when

0 < prevalence < 1
results conditionally independent given the status

Example

m.ask("Repeated testing: prevalence 0.01, sensitivity 0.9, specificity 0.95, positives = 2, negatives = 0. Posterior probability of the disease?")

→ posterior = 0.765957, posterior_odds = 3.27273

Posterior over several discrete hypothesesdiscrete_bayes_posterior

Posterior probabilities of k hypotheses from their priors and the likelihood of the observed data under each; the evidence and the most probable hypothesis (MAP).

P(H_i\mid D) = \frac{P(D\mid H_i)\,P(H_i)}{\sum_j P(D\mid H_j)\,P(H_j)}

Inputs

priors prior probabilities of the hypotheses (sum to 1) vector
likelihoods likelihood of the data under each hypothesis vector

Returns

posterior posterior probabilities
evidence marginal probability of the data
map_hypothesis number (1-based) of the most probable hypothesis
map_probability its posterior probability

Applies when

priors sum to 1
likelihoods >= 0, not all zero where the prior is positive

Example

m.ask("Posterior probabilities with priors = [0.5, 0.3, 0.2] and likelihoods = [0.1, 0.4, 0.7].")

→ posterior = [0.161290322581, 0.387096774194, 0.451612903226], evidence = 0.31, map_hypothesis = 3, map_probability = 0.451613

Conjugate priors 17

Bayesian A/B test: probability that A beats Bab_test_probability_beta

P(p_A > p_B) for independent Beta posteriors of two success probabilities (alpha_b a whole number, as after a uniform or any integer-parameter prior).

P(p_B > p_A) = \sum_{i=0}^{\alpha_B - 1}\frac{B(\alpha_A + i,\ \beta_A + \beta_B)}{(\beta_B + i)\,B(1 + i,\ \beta_B)\,B(\alpha_A, \beta_A)}

Inputs

alpha_a posterior alpha of A
beta_a posterior beta of A
alpha_b posterior alpha of B (whole number)
beta_b posterior beta of B

Returns

prob_a_greater P(p_A > p_B)
prob_b_greater P(p_B > p_A)

Applies when

all parameters > 0
alpha_b a whole number

Example

m.ask("Bayesian A/B test: alpha_a = 1, beta_a = 1, alpha_b = 1, beta_b = 1. Probability that A beats B?")

→ prob_a_greater = 0.5, prob_b_greater = 0.5

Beta-binomial conjugate posteriorbeta_binomial_posterior

Update a Beta(alpha, beta) prior on a success probability with k successes in n trials: posterior parameters, mean, variance and standard deviation.

p\mid k \sim \mathrm{Beta}(\alpha + k,\ \beta + n - k),\quad E[p\mid k] = \frac{\alpha + k}{\alpha + \beta + n}

Inputs

alpha prior Beta alpha (> 0)
beta prior Beta beta (> 0)
successes number of successes k integer
trials number of trials n integer

Returns

post_alpha posterior alpha = alpha + k
post_beta posterior beta = beta + n - k
posterior_mean posterior mean
posterior_variance posterior variance
posterior_sd posterior standard deviation

Applies when

alpha, beta > 0
0 <= k <= n

Example

m.ask("Beta-binomial posterior: prior Beta(2, 2), 7 successes in 10 trials. Find the posterior mean.")

→ post_alpha = 9, post_beta = 5, posterior_mean = 0.642857, posterior_variance = 0.0153061, posterior_sd = 0.123718

Dirichlet-multinomial conjugate posteriordirichlet_multinomial_posterior

Update a Dirichlet(alpha_1, ..., alpha_k) prior on category probabilities with category counts: posterior concentrations, means and variances.

\theta\mid n \sim \mathrm{Dir}(\alpha_1 + n_1, \ldots, \alpha_k + n_k),\quad E[\theta_j] = \frac{\alpha_j + n_j}{\alpha_0 + n}

Inputs

alpha prior concentrations alpha_j (> 0) vector
counts observed category counts vector

Returns

post_alpha alpha_j + n_j
posterior_means posterior means (alpha_j + n_j)/(alpha_0 + n)
posterior_variances posterior variances

Applies when

alpha_j > 0
counts whole numbers >= 0
same length

Example

m.ask("Dirichlet-multinomial posterior with alpha = [1, 1, 1] and counts = [10, 5, 3]: posterior means?")

→ post_alpha = [11.0, 6.0, 4.0], posterior_means = [0.52380952381, 0.285714285714, 0.190476190476], posterior_variances = [0.0113378684807, 0.00927643784787, 0.00700886415172]

Dirichlet-multinomial probability of a count vector (evidence)dirichlet_multinomial_evidence

Prior predictive probability (marginal likelihood) of observed category counts when the category probabilities have a Dirichlet prior.

P(n\mid\alpha) = \frac{n!\,\Gamma(\alpha_0)}{\Gamma(n + \alpha_0)}\prod_j\frac{\Gamma(n_j + \alpha_j)}{n_j!\,\Gamma(\alpha_j)}

Inputs

alpha prior concentrations (> 0) vector
counts category counts vector

Returns

probability P(counts | alpha)
log_probability its natural log

Applies when

alpha_j > 0
counts whole numbers >= 0

Example

m.ask("Dirichlet-multinomial evidence for counts = [2, 1, 0] with alpha = [1, 1, 1].")

→ probability = 0.1, log_probability = -2.30259

Gamma-exponential conjugate posterior of a rategamma_exponential_posterior

Update a Gamma(shape, rate) prior on the rate of exponential waiting times with observed times.

\lambda\mid t \sim \mathrm{Gamma}\!\left(\alpha + n,\ \beta + \textstyle\sum t_i\right)

Inputs

shape prior shape (> 0)
rate prior rate (> 0)
times observed waiting times (> 0) vector

Returns

post_shape shape + n
post_rate rate + sum of times
posterior_mean posterior mean of the rate
posterior_variance posterior variance
posterior_sd posterior standard deviation

Applies when

shape, rate > 0
times > 0

Example

m.ask("Gamma-exponential posterior: prior Gamma(2, 1), exponential waiting times = [0.5, 1.2, 0.8]. Posterior mean of the rate?")

→ post_shape = 5, post_rate = 3.5, posterior_mean = 1.42857, posterior_variance = 0.408163, posterior_sd = 0.638877

Gamma-Poisson conjugate posterior of a rategamma_poisson_posterior

Update a Gamma(shape, rate) prior on a Poisson rate with observed counts: posterior parameters, mean, variance and standard deviation.

\lambda\mid x \sim \mathrm{Gamma}\!\left(\alpha + \textstyle\sum x_i,\ \beta + n\right)

Inputs

shape prior shape alpha (> 0)
rate prior rate beta (> 0)
counts observed counts (whole numbers), one per unit interval vector

Returns

post_shape alpha + sum of counts
post_rate beta + number of intervals
posterior_mean posterior mean of the rate
posterior_variance posterior variance
posterior_sd posterior standard deviation

Applies when

shape, rate > 0
counts whole numbers >= 0

Example

m.ask("Gamma-Poisson posterior: prior Gamma(2, 1), counts = [3, 5, 4]. Posterior mean of the rate?")

→ post_shape = 14, post_rate = 4, posterior_mean = 3.5, posterior_variance = 0.875, posterior_sd = 0.935414

Laplace's rule of successionrule_of_succession

Probability that the next trial is a success after k successes in n trials, with a uniform prior on the success probability: (k + 1)/(n + 2).

P(\text{next success}) = \frac{k+1}{n+2}

Inputs

successes number of successes k integer
trials number of trials n integer

Returns

p_next posterior predictive probability that the next trial succeeds
mle observed proportion k/n (none when n = 0)

Applies when

0 <= k <= n
uniform Beta(1, 1) prior

Example

m.ask("Rule of succession: 7 successes in 10 trials. Probability the next trial is a success?")

→ p_next = 0.666667, mle = 0.7

MAP estimate of a Poisson rate (gamma prior)map_gamma_poisson

Posterior mode of a Poisson rate with a Gamma(shape, rate) prior and observed counts, with the MLE and the posterior mean for comparison.

\hat\lambda_{MAP} = \frac{\alpha + \sum x_i - 1}{\beta + n}

Inputs

shape prior shape (> 0)
rate prior rate (> 0)
counts observed counts vector

Returns

map (shape + sum - 1)/(rate + n)
mle sample mean of the counts
posterior_mean (shape + sum)/(rate + n)

Applies when

shape + sum of counts > 1

Example

m.ask("MAP estimate of the Poisson rate with prior Gamma(2, 1) and counts = [3, 5, 4].")

→ map = 3.25, mle = 4, posterior_mean = 3.5

MAP estimate of a proportion (beta-binomial)map_beta_binomial

The posterior mode (maximum a posteriori estimate) of p after k successes in n trials with a Beta(alpha, beta) prior, with the maximum-likelihood estimate and the posterior mean for comparison.

\hat p_{MAP} = \frac{\alpha + k - 1}{\alpha + \beta + n - 2}

Inputs

alpha prior Beta alpha (> 0)
beta prior Beta beta (> 0)
successes number of successes k integer
trials number of trials n integer

Returns

map posterior mode (alpha + k - 1)/(alpha + beta + n - 2)
mle k/n
posterior_mean (alpha + k)/(alpha + beta + n)

Applies when

alpha + k > 1 and beta + n - k > 1 (interior mode)

Example

m.ask("MAP estimate of p with prior Beta(2, 2) and 7 successes in 10 trials.")

→ map = 0.666667, mle = 0.7, posterior_mean = 0.642857

Normal mean with known variance: posterior and credible interval from a sample meannormal_normal_summary

Posterior mean, standard deviation, credible interval and shrinkage factor of a normal mean from the sample mean of n observations with known sigma and a normal prior.

\mu_n = B\mu_0 + (1-B)\bar x,\ B = \frac{\sigma^2/n}{\sigma^2/n + \tau_0^2},\quad \mu_n \pm z_{(1+c)/2}\,\tau_n

Inputs

prior_mean prior mean
prior_sd prior standard deviation (> 0)
sigma known standard deviation of one observation (> 0)
sample_mean sample mean
n sample size (>= 1) integer
level credibility level as a fraction (0.95 for 95%; default 0.95) optional

Returns

posterior_mean posterior mean
posterior_sd posterior standard deviation
lower lower credible limit
upper upper credible limit
shrinkage weight on the prior mean (sigma^2/n)/(sigma^2/n + tau0^2)

Applies when

prior sd > 0
sigma > 0
n >= 1
0 < level < 1

Example

m.ask("Normal prior with prior mean = 0, prior sd = 1; sample mean = 1 from n = 4 with known sigma = 2. Posterior mean and 95% credible interval?")

→ posterior_mean = 0.5, posterior_sd = 0.707107, lower = -0.885904, upper = 1.8859, shrinkage = 0.5

Normal mean with known variance: posterior from datanormal_normal_posterior

Posterior of a normal mean theta with a N(prior mean, prior sd^2) prior, given observations with known standard deviation sigma.

\frac{1}{\tau_n^2} = \frac{1}{\tau_0^2} + \frac{n}{\sigma^2},\quad \mu_n = \tau_n^2\left(\frac{\mu_0}{\tau_0^2} + \frac{n\bar x}{\sigma^2}\right)

Inputs

prior_mean prior mean mu0
prior_sd prior standard deviation tau0 (> 0)
sigma known standard deviation of each observation (> 0)
data the observations vector

Returns

posterior_mean posterior mean
posterior_sd posterior standard deviation
posterior_variance posterior variance

Applies when

prior sd > 0
sigma > 0 known

Example

m.ask("Normal-normal posterior: prior mean = 10, prior sd = 2, known sigma = 3, data = [12, 9, 14]. Posterior mean?")

→ posterior_mean = 10.9524, posterior_sd = 1.30931, posterior_variance = 1.71429

Normal model with unknown mean and variance (normal-inverse-gamma)normal_inverse_gamma_posterior

Conjugate normal-inverse-gamma update for normal data with unknown mean and variance: posterior parameters, posterior mean of the variance and the Student-t credible interval for the mean.

\beta_n = \beta_0 + \tfrac12\sum(x_i-\bar x)^2 + \frac{\kappa_0 n(\bar x-\mu_0)^2}{2(\kappa_0+n)},\quad \mu\mid x \sim t_{2\alpha_n}\!\left(\mu_n,\ \frac{\beta_n}{\alpha_n\kappa_n}\right)

Inputs

mu0 prior location of the mean
kappa0 prior pseudo-observations for the mean (> 0)
alpha0 prior inverse-gamma shape (> 0)
beta0 prior inverse-gamma scale (> 0)
data the observations vector
level credibility level as a fraction (0.95 for 95%; default 0.95) optional

Returns

mu_n posterior location
kappa_n kappa0 + n
alpha_n alpha0 + n/2
beta_n posterior scale
sigma2_mean posterior mean of sigma^2 (alpha_n > 1)
lower lower credible limit for the mean
upper upper credible limit for the mean

Applies when

kappa0, alpha0, beta0 > 0
0 < level < 1

Example

m.ask("Normal-inverse-gamma posterior: mu0 = 0, kappa0 = 1, alpha0 = 2, beta0 = 2, data = [1, -1]. 95% credible interval for the mean?")

→ mu_n = 0, kappa_n = 3, alpha_n = 3, beta_n = 3, sigma2_mean = 1.5, lower = -1.41273, upper = 1.41273

Posterior predictive count (negative binomial)gamma_poisson_predictive

Probability of exactly k events in the next unit interval when the Poisson rate has a Gamma(shape, rate) posterior; the cumulative probability and the predictive mean and variance.

P(X = k) = \frac{\Gamma(\alpha + k)}{\Gamma(\alpha)\,k!}\left(\frac{\beta}{\beta+1}\right)^{\alpha}\left(\frac{1}{\beta+1}\right)^{k}

Inputs

shape posterior shape (> 0)
rate posterior rate (> 0)
events number of events k in the next interval integer

Returns

probability P(exactly k events)
prob_at_most P(at most k events)
predictive_mean shape/rate
predictive_variance shape/rate (1 + 1/rate)

Applies when

shape, rate > 0
k whole number >= 0

Example

m.ask("Posterior Gamma(1, 1) for a Poisson rate: probability of exactly 1 events in the next interval (posterior predictive)?")

→ probability = 0.25, prob_at_most = 0.75, predictive_mean = 1, predictive_variance = 2

Posterior predictive distribution of a new observationnormal_posterior_predictive

Mean, standard deviation and tail probabilities of the next observation in the normal model with known sigma, after a sample mean of n observations and a normal prior.

Y_{new}\mid \bar x \sim N\!\left(\mu_n,\ \tau_n^2 + \sigma^2\right)

Inputs

prior_mean prior mean
prior_sd prior standard deviation (> 0)
sigma known standard deviation (> 0)
sample_mean sample mean
n sample size (>= 1) integer
threshold value c for P(Y > c)

Returns

predictive_mean posterior mean
predictive_sd sqrt(posterior variance + sigma^2)
prob_above P(Y_new > c)
prob_below P(Y_new < c)

Applies when

prior sd > 0
sigma > 0
n >= 1

Example

m.ask("Posterior predictive: prior mean = 0, prior sd = 1, known sigma = 2, sample mean = 1, n = 4. Probability that the next observation exceeds 0.5?")

→ predictive_mean = 0.5, predictive_sd = 2.12132, prob_above = 0.5, prob_below = 0.5

Posterior predictive probability (beta-binomial)beta_binomial_predictive

Probability of exactly j successes in m future trials when the success probability has a Beta(alpha, beta) posterior; the cumulative probability and the predictive mean.

P(J = j) = \binom{m}{j}\frac{B(\alpha + j,\ \beta + m - j)}{B(\alpha, \beta)}

Inputs

alpha posterior alpha (> 0)
beta posterior beta (> 0)
future_trials number of future trials m (>= 1) integer
future_successes number of future successes j integer

Returns

probability P(exactly j successes in m trials)
prob_at_most P(at most j)
prob_more P(more than j)
predictive_mean expected number of successes m alpha/(alpha+beta)

Applies when

alpha, beta > 0
0 <= j <= m

Example

m.ask("Posterior Beta(1, 1): probability of exactly 2 successes in the next 4 trials (posterior predictive)?")

→ probability = 0.2, prob_at_most = 0.6, prob_more = 0.4, predictive_mean = 2

Posterior probability that p exceeds a valuebeta_posterior_tail_probability

P(p > c) and P(p < c) under a Beta(alpha, beta) posterior, optionally after k successes in n trials.

P(p > c \mid k) = 1 - I_c(\alpha + k,\ \beta + n - k)

Inputs

alpha Beta alpha (> 0)
beta Beta beta (> 0)
threshold the value c, 0 < c < 1
successes successes to add (default 0) integer optional
trials trials to add (default 0) integer optional

Returns

prob_above posterior probability P(p > c)
prob_below posterior probability P(p < c)

Applies when

alpha, beta > 0
0 < c < 1
0 <= k <= n

Example

m.ask("For a Beta(2, 1) posterior, what is the probability that p exceeds 0.5?")

→ prob_above = 0.75, prob_below = 0.25

Probability that one Poisson rate exceeds another (gamma posteriors)probability_rate_greater_gamma

P(lambda_1 > lambda_2) for independent Gamma(shape, rate) posteriors of two Poisson rates, and the ratio of their posterior means.

P(\lambda_1 > \lambda_2) = 1 - I_{b_1/(b_1+b_2)}(a_1, a_2)

Inputs

shape1 posterior shape of rate 1
rate1 posterior rate parameter of rate 1
shape2 posterior shape of rate 2
rate2 posterior rate parameter of rate 2

Returns

prob_1_greater P(lambda_1 > lambda_2)
prob_2_greater P(lambda_2 > lambda_1)
rate_ratio_mean ratio of posterior means

Applies when

all parameters > 0

Example

m.ask("Compare gamma posteriors shape1 = 5, rate1 = 1, shape2 = 5, rate2 = 1.")

→ prob_1_greater = 0.5, prob_2_greater = 0.5, rate_ratio_mean = 1

Credible intervals 7

Credible interval for a Poisson rategamma_poisson_credible_interval

Equal-tailed credible interval of a Poisson rate after observed counts with a Gamma(shape, rate) prior.

\left[\frac{G^{-1}_{\alpha'}\!\left(\frac{1-c}{2}\right)}{\beta'},\ \frac{G^{-1}_{\alpha'}\!\left(\frac{1+c}{2}\right)}{\beta'}\right]

Inputs

shape prior shape (> 0)
rate prior rate (> 0)
counts observed counts vector
level credibility level as a fraction (0.95 for 95%; default 0.95) optional

Returns

post_shape posterior shape
post_rate posterior rate
lower lower limit
upper upper limit

Applies when

shape, rate > 0
counts whole numbers >= 0
0 < level < 1

Example

m.ask("90% credible interval for a Poisson rate: prior Gamma(1, 1), counts = [0].")

→ post_shape = 1, post_rate = 2, lower = 0.0256466, upper = 1.49787

Credible interval for a proportion (beta-binomial)beta_binomial_credible_interval

Equal-tailed credible interval of a success probability after k successes in n trials with a Beta(alpha, beta) prior.

p\mid k \sim \mathrm{Beta}(\alpha+k,\beta+n-k),\quad [q_{(1-c)/2},\ q_{(1+c)/2}]

Inputs

alpha prior Beta alpha (> 0)
beta prior Beta beta (> 0)
successes number of successes k integer
trials number of trials n integer
level credibility level as a fraction (0.95 for 95%; default 0.95) optional

Returns

post_alpha posterior alpha
post_beta posterior beta
lower lower credible limit
upper upper credible limit

Applies when

alpha, beta > 0
0 <= k <= n
0 < level < 1

Example

m.ask("90% credible interval for p: prior Beta(1, 1), 1 successes in 1 trials.")

→ post_alpha = 2, post_beta = 1, lower = 0.223607, upper = 0.974679

Equal-tailed credible interval of a beta distributionbeta_credible_interval

The central (equal-tailed) credible interval and the median of a Beta(alpha, beta) posterior.

\left[F^{-1}_{\mathrm{Beta}(\alpha,\beta)}\!\left(\tfrac{1-c}{2}\right),\ F^{-1}_{\mathrm{Beta}(\alpha,\beta)}\!\left(\tfrac{1+c}{2}\right)\right]

Inputs

alpha Beta alpha (> 0)
beta Beta beta (> 0)
level credibility level as a fraction (0.95 for 95%; default 0.95) optional

Returns

lower lower limit (quantile (1 - level)/2)
upper upper limit (quantile (1 + level)/2)
median posterior median

Applies when

alpha, beta > 0
0 < level < 1

Example

m.ask("Equal-tailed 90% credible interval for a Beta(1, 1) posterior.")

→ lower = 0.05, upper = 0.95, median = 0.5

Equal-tailed credible interval of a gamma distributiongamma_credible_interval

The central credible interval and median of a Gamma(shape, rate) posterior.

\left[G^{-1}\!\left(\tfrac{1-c}{2}\right),\ G^{-1}\!\left(\tfrac{1+c}{2}\right)\right]

Inputs

shape shape (> 0)
rate rate (> 0)
level credibility level as a fraction (0.95 for 95%; default 0.95) optional

Returns

lower lower limit
upper upper limit
median median

Applies when

shape, rate > 0
0 < level < 1

Example

m.ask("90% credible interval of a Gamma(1, 2) posterior.")

→ lower = 0.0256466, upper = 1.49787, median = 0.346574

Highest-density interval of a beta posteriorbeta_hdi

The shortest interval holding a stated posterior probability of a unimodal Beta(alpha, beta): equal density at both ends.

F(u) - F(l) = c,\quad f(l) = f(u)

Inputs

alpha Beta alpha (> 1)
beta Beta beta (> 1)
level credibility level as a fraction (0.95 for 95%; default 0.95) optional

Returns

lower lower limit
upper upper limit
width interval width

Applies when

alpha > 1 and beta > 1 (unimodal, zero density at both ends)
0 < level < 1

Example

m.ask("95% highest-density interval of a Beta(3, 9) posterior.")

→ lower = 0.0405552, upper = 0.483724, width = 0.443168

Highest-density interval of a gamma posteriorgamma_hdi

The shortest interval holding a stated probability of a unimodal Gamma(shape, rate): equal density at both ends.

G(u) - G(l) = c,\quad g(l) = g(u)

Inputs

shape shape (> 1)
rate rate (> 0)
level credibility level as a fraction (0.95 for 95%; default 0.95) optional

Returns

lower lower limit
upper upper limit
width interval width

Applies when

shape > 1
rate > 0
0 < level < 1

Example

m.ask("95% highest-density interval of a Gamma(5, 2) posterior.")

→ lower = 0.60348, upper = 4.71511, width = 4.11163

Jeffreys interval for a binomial proportionjeffreys_interval

The equal-tailed Bayesian interval for a proportion under the Jeffreys Beta(1/2, 1/2) prior (lower limit 0 when k = 0 and upper limit 1 when k = n).

\left[B^{-1}\!\left(\tfrac{1-c}{2};\,k+\tfrac12,\,n-k+\tfrac12\right),\ B^{-1}\!\left(\tfrac{1+c}{2};\,k+\tfrac12,\,n-k+\tfrac12\right)\right]

Inputs

successes number of successes k integer
trials number of trials n (>= 1) integer
level credibility level as a fraction (0.95 for 95%; default 0.95) optional

Returns

lower lower limit
upper upper limit
posterior_mean (k + 1/2)/(n + 1)

Applies when

0 <= k <= n
n >= 1
0 < level < 1

Example

m.ask("Jeffreys interval at 95% for 7 successes in 50 trials.")

→ lower = 0.0648717, upper = 0.255325, posterior_mean = 0.147059

Decision analysis 9

Bayes action: minimum posterior expected lossbayes_action_minimum_expected_loss

Posterior expected loss of each action from a loss matrix and (posterior) state probabilities; the Bayes action minimizes it.

\rho(a) = \sum_s P(s\mid x)\,L(a,s),\quad a^* = \arg\min_a \rho(a)

Inputs

losses loss matrix: one row per action, one column per state matrix
probabilities posterior state probabilities (sum to 1) vector

Returns

expected_losses expected loss of each action
bayes_action action with least expected loss
bayes_risk its expected loss

Applies when

probabilities sum to 1

Example

m.ask("Bayes action: losses = [[0, 10], [5, 0]], probabilities = [0.3, 0.7]. Minimum expected loss?")

→ expected_losses = [7.0, 1.5], bayes_action = 2, bayes_risk = 1.5

Best action by expected monetary valueexpected_value_decision

Expected payoff of each action over uncertain states with known probabilities, and the action with the highest expected value.

EV(a) = \sum_s P(s)\,u(a, s),\quad a^* = \arg\max_a EV(a)

Inputs

payoffs payoff matrix: one row per action, one column per state matrix
probabilities state probabilities (sum to 1) vector

Returns

expected_values expected payoff of each action
best_action number (1-based) of the best action (first if tied)
best_expected_value its expected payoff

Applies when

probabilities sum to 1
one probability per column

Example

m.ask("Expected value decision: payoffs = [[100, -20], [40, 30]], probabilities = [0.4, 0.6]. Best action?")

→ expected_values = [28.0, 34.0], best_action = 2, best_expected_value = 34

Break-even (indifference) probability of two actionsbreak_even_probability

The probability of state 1 at which two actions have equal expected value, and which action is preferred above it.

p^* = \frac{u_{22} - u_{12}}{(u_{11} - u_{12}) - (u_{21} - u_{22})}

Inputs

payoffs 2 x 2 payoff matrix [[action 1 in state 1, in state 2], [action 2 ...]] matrix

Returns

break_even probability of state 1 at which the expected values are equal
preferred_above action preferred when P(state 1) exceeds it
expected_value_at_break_even the common expected value there

Applies when

a crossing strictly between 0 and 1

Example

m.ask("Break-even probability for payoffs = [[100, -20], [40, 30]].")

→ break_even = 0.454545, preferred_above = 1, expected_value_at_break_even = 34.5455

Certainty equivalent with logarithmic utilitycertainty_equivalent_log_utility

Expected log utility, certainty equivalent (the weighted geometric mean) and risk premium of a lottery over positive wealth levels.

CE = \exp\!\left(\sum_i p_i \ln x_i\right) = \prod_i x_i^{p_i}

Inputs

outcomes positive wealth levels vector
probabilities their probabilities (sum to 1) vector

Returns

expected_log_utility E[ln x]
certainty_equivalent exp(E[ln x])
risk_premium E[x] - CE
expected_value E[x]

Applies when

outcomes > 0
probabilities sum to 1

Example

m.ask("Logarithmic utility: outcomes = [100, 400], probabilities = [0.5, 0.5]. Certainty equivalent?")

→ expected_log_utility = 5.29832, certainty_equivalent = 200, risk_premium = 50, expected_value = 250

Decision threshold probability (treat / do not treat)treatment_threshold

The probability above which acting (treating) has lower expected loss than not acting, from the loss of a false positive and of a false negative; the decision for a given probability.

p^* = \frac{L_{FP}}{L_{FP} + L_{FN}}

Inputs

loss_fp loss of acting when the condition is absent (false positive)
loss_fn loss of not acting when the condition is present (false negative)
probability probability of the condition (optional) optional

Returns

threshold loss_fp / (loss_fp + loss_fn)
decision treat / do not treat (if p given)

Applies when

losses > 0

Example

m.ask("Treatment threshold with loss_fp = 1 and loss_fn = 9; the probability of disease is 0.05. Treat?")

→ threshold = 0.1, decision = do not treat

Decision under uncertainty: maximin, maximax, Laplace, minimax regret, Hurwiczdecision_criteria_under_uncertainty

The best action and its criterion value when state probabilities are unknown, by a stated criterion.

\text{maximin: }\max_a\min_s u,\quad \text{regret } r(a,s) = \max_{a'}u(a',s) - u(a,s)

Inputs

payoffs payoff matrix: one row per action, one column per state matrix
criterion maximin, maximax, laplace, minimax regret or hurwicz choice
optimism Hurwicz optimism coefficient alpha (Hurwicz only) optional

Returns

scores criterion score of each action (maximum regret for minimax regret)
best_action number of the chosen action (first if tied)
criterion_value criterion value of the chosen action

Applies when

0 <= alpha <= 1 for Hurwicz

Example

m.ask("Use the maximin criterion with payoffs = [[100, -20], [40, 30], [60, 10]].")

→ scores = [-20.0, 30.0, 10.0], best_action = 2, criterion_value = 30

Expected utility and certainty equivalent (exponential utility)certainty_equivalent_exponential_utility

Expected utility, certainty equivalent and risk premium of a lottery for a risk-averse decision maker with exponential utility u(x) = 1 - exp(-x/R).

CE = -R\ln\sum_i p_i e^{-x_i/R}

Inputs

outcomes monetary outcomes vector
probabilities their probabilities (sum to 1) vector
risk_tolerance risk tolerance R (> 0)

Returns

expected_utility E[1 - exp(-x/R)]
certainty_equivalent -R ln E[exp(-x/R)]
risk_premium expected value - certainty equivalent
expected_value expected value

Applies when

probabilities sum to 1
R > 0

Example

m.ask("Exponential utility with risk tolerance = 500: outcomes = [1000, 0], probabilities = [0.5, 0.5]. Certainty equivalent?")

→ expected_utility = 0.432332, certainty_equivalent = 283.11, risk_premium = 216.89, expected_value = 500

Expected value of perfect information (EVPI)expected_value_of_perfect_information

The most a decision maker should pay to learn the true state before acting: expected value with perfect information minus the best expected value without it.

EVPI = \sum_s P(s)\max_a u(a,s) - \max_a\sum_s P(s)\,u(a,s)

Inputs

payoffs payoff matrix: one row per action, one column per state matrix
probabilities state probabilities (sum to 1) vector

Returns

evpi expected value of perfect information
ev_with_perfect_information expected value of the best action in each state
best_expected_value best expected value without information

Applies when

probabilities sum to 1

Example

m.ask("Expected value of perfect information: payoffs = [[100, -20], [40, 30]], probabilities = [0.4, 0.6].")

→ evpi = 24, ev_with_perfect_information = 58, best_expected_value = 34

Expected value of sample information (EVSI)expected_value_of_sample_information

Value of an imperfect test before acting: for each possible signal, update the state probabilities by Bayes' theorem, choose the best action, and compare the expected value with that of acting now.

EVSI = \sum_z P(z)\max_a\sum_s P(s\mid z)\,u(a,s) - \max_a\sum_s P(s)\,u(a,s)

Inputs

payoffs payoff matrix: one row per action, one column per state matrix
priors prior state probabilities (sum to 1) vector
likelihoods likelihood matrix P(signal | state): one row per signal, one column per state matrix

Returns

evsi expected value of sample information
ev_with_information expected value with the test
ev_without_information best expected value without the test
signal_probabilities probability of each signal
best_action_per_signal best action after each signal

Applies when

priors sum to 1
each likelihood column sums to 1

Example

m.ask("Find the EVSI with payoff matrix [[50, -10], [0, 0]], prior probabilities [0.5, 0.5] and likelihood matrix [[0.9, 0.2], [0.1, 0.8]].")

→ evsi = 1.5, ev_with_information = 21.5, ev_without_information = 20, signal_probabilities = [0.55, 0.45], best_action_per_signal = [1, 2]

Information and quadrature for Bayesian computation 3

Gauss-Hermite rule for expectations under a normal distributiongauss_hermite_rule

Nodes and weights of the n-point Gauss-Hermite rule for E[f(Z)], Z standard normal (probabilists' Hermite polynomials; weights sum to 1; exact for polynomials of degree up to 2n - 1). For N(mu, sigma^2) use mu + sigma x_i.

E[f(Z)] \approx \sum_{i=1}^{n} w_i f(x_i),\quad He_n(x_i) = 0,\quad w_i = \frac{n!}{n^2\,He_{n-1}(x_i)^2}

Inputs

n number of nodes (1 to 40) integer

Returns

nodes nodes x_i (ascending)
weights weights w_i (sum to 1)

Applies when

1 <= n <= 40

Example

m.ask("Gauss-Hermite rule with n = 3 nodes for a normal expectation.")

→ nodes = [-1.73205080757, 0.0, 1.73205080757], weights = [0.166666666667, 0.666666666667, 0.166666666667]

Kullback-Leibler divergence between two beta distributionskl_divergence_beta

KL(Beta(alpha1, beta1) || Beta(alpha2, beta2)): for example the information gained from a beta prior (distribution 2) to a beta posterior (distribution 1).

KL = \ln\frac{B(\alpha_2,\beta_2)}{B(\alpha_1,\beta_1)} + (\alpha_1-\alpha_2)\psi(\alpha_1) + (\beta_1-\beta_2)\psi(\beta_1) + (\alpha_2-\alpha_1+\beta_2-\beta_1)\psi(\alpha_1+\beta_1)

Inputs

alpha1 alpha of distribution 1
beta1 beta of distribution 1
alpha2 alpha of distribution 2
beta2 beta of distribution 2

Returns

kl KL divergence in nats

Applies when

all parameters > 0

Example

m.ask("KL divergence between beta distributions: alpha1 = 2, beta1 = 1, alpha2 = 1, beta2 = 1.")

→ kl = 0.193147

Kullback-Leibler divergence between two normal distributionskl_divergence_normal

KL(N1 || N2) and KL(N2 || N1), e.g. the information gained from a normal prior to a normal posterior.

KL(N_1\Vert N_2) = \ln\frac{\sigma_2}{\sigma_1} + \frac{\sigma_1^2 + (\mu_1-\mu_2)^2}{2\sigma_2^2} - \frac12

Inputs

mean1 mean of distribution 1
sd1 sd of distribution 1 (> 0)
mean2 mean of distribution 2
sd2 sd of distribution 2 (> 0)

Returns

kl_1_to_2 KL(N1 || N2) in nats
kl_2_to_1 KL(N2 || N1) in nats

Applies when

sd1, sd2 > 0

Example

m.ask("KL divergence between normal distributions: mean1 = 0, sd1 = 1, mean2 = 1, sd2 = 2.")

→ kl_1_to_2 = 0.443147, kl_2_to_1 = 1.30685

Monte Carlo estimation 6

Control-variate Monte Carlo estimatecontrol_variate_estimate

Variance-reduced estimate of E[Y] from paired draws (y_i, c_i) of the output and a control variate with known mean: optimal coefficient, estimate and the fraction of variance left (1 - r^2).

\hat\mu_{cv} = \bar y - \hat b\,(\bar c - \mu_c),\quad \hat b = \frac{\sum (y_i-\bar y)(c_i-\bar c)}{\sum (c_i-\bar c)^2}

Inputs

y simulated outputs y_i vector
c control values c_i (same draws) vector
control_mean known mean of the control

Returns

estimate mean(y) - b (mean(c) - mu_c)
coefficient b = cov(y, c)/var(c)
plain_estimate mean(y)
variance_reduction variance factor 1 - r^2

Applies when

outputs and control vary
n >= 3

Example

m.ask("Control variate: outputs = [1, 2, 3], controls = [1, 2, 3], known mean = 2.5.")

→ estimate = 2.5, coefficient = 1, plain_estimate = 2, variance_reduction = 0

Importance sampling estimate and effective sample sizeimportance_sampling_estimate

Self-normalized and plain importance-sampling estimates of E_p[f] from draws of a proposal q, with the effective sample size of the weights p/q.

\hat\mu = \frac{\sum_i w_i f(x_i)}{\sum_i w_i},\ w_i = \frac{p(x_i)}{q(x_i)},\quad ESS = \frac{(\sum w_i)^2}{\sum w_i^2}

Inputs

values f(x_i) at the draws vector
target_density target density p(x_i) (or unnormalized) vector
proposal_density proposal density q(x_i) > 0 vector

Returns

estimate self-normalized estimate sum w f / sum w
unnormalized_estimate mean of w f (needs normalized densities)
effective_sample_size (sum w)^2 / sum w^2

Applies when

q > 0
p >= 0, not all zero

Example

m.ask("Importance sampling estimate with f values [0, 1], target density [0.5, 0.5] and proposal density [0.25, 0.75].")

→ estimate = 0.25, unnormalized_estimate = 0.333333, effective_sample_size = 1.6

Monte Carlo estimate of an integral from uniform samplesmonte_carlo_integration

Estimate of the integral of f over [a, b] from f evaluated at n points drawn uniformly on [a, b], with its Monte Carlo standard error.

\hat I = \frac{b-a}{n}\sum_i f(x_i),\quad SE = (b-a)\frac{s_f}{\sqrt n}

Inputs

values f(x_i) at uniform sample points vector
lower lower limit a
upper upper limit b (> a)

Returns

estimate (b - a) times the mean of f
standard_error (b - a) s / sqrt(n)
n number of samples

Applies when

b > a
n >= 2

Example

m.ask("Monte Carlo integration with values = [1, 3]: lower = 0, upper = 2.")

→ estimate = 4, standard_error = 2, n = 2

Number of Monte Carlo draws for a target accuracymonte_carlo_sample_size

How many independent draws make the half-width of the confidence interval of a Monte Carlo mean at most the stated value, given the standard deviation of one draw.

n = \left\lceil\left(\frac{z_{(1+c)/2}\,\sigma}{E}\right)^2\right\rceil

Inputs

sigma standard deviation of one draw (> 0)
half_width target half-width (> 0)
level credibility level as a fraction (0.95 for 95%; default 0.95) optional

Returns

draws required number of draws (rounded up)
z normal quantile z_{(1+c)/2}
standard_error_target target standard error half_width/z

Applies when

sigma, half-width > 0
0 < level < 1

Example

m.ask("Monte Carlo sample size: standard deviation 10, half-width 0.5, 95% confidence. How many draws?")

→ draws = 1537, z = 1.95996, standard_error_target = 0.255107

Point of the Halton low-discrepancy sequencehalton_point

Coordinates of the point with a given index in the (unscrambled) Halton sequence of a given dimension: the radical inverses of the index in the first prime bases (index 0 is the origin).

\phi_b(i) = \sum_k d_k\,b^{-k-1}\ \text{ for } i = \sum_k d_k b^k

Inputs

index index of the point (0, 1, 2, ...) integer
dimension dimension (1 to 20) integer

Returns

point coordinates in [0, 1)
bases prime bases used

Applies when

index >= 0
1 <= dimension <= 20

Example

m.ask("Halton sequence point with index = 5 and dimension = 3.")

→ point = [0.625, 0.777777777778, 0.04], bases = [2.0, 3.0, 5.0]

Stratified sampling estimate and standard errorstratified_sampling_estimate

Estimate of a mean from strata with known weights, the sample mean, standard deviation and size of each stratum, and its standard error.

\hat\mu = \sum_h W_h\bar y_h,\quad SE = \sqrt{\sum_h W_h^2 s_h^2/n_h}

Inputs

weights stratum weights (sum to 1) vector
means stratum sample means vector
sds stratum sample standard deviations vector
sizes stratum sample sizes vector

Returns

estimate sum W_h mean_h
standard_error sqrt(sum W_h^2 s_h^2 / n_h)

Applies when

weights sum to 1
sizes >= 1

Example

m.ask("Stratified sampling: weights = [0.6, 0.4], means = [10, 20], sds = [2, 5], sizes = [30, 20].")

→ estimate = 14, standard_error = 0.497996

Odds and likelihood ratios 4

Likelihood ratios of a diagnostic testlikelihood_ratios

Positive and negative likelihood ratios and the diagnostic odds ratio from sensitivity and specificity.

LR^{+} = \frac{\mathrm{Se}}{1-\mathrm{Sp}},\quad LR^{-} = \frac{1-\mathrm{Se}}{\mathrm{Sp}}

Inputs

sensitivity sensitivity P(positive | condition)
specificity specificity P(negative | no condition)

Returns

lr_positive LR+ = sensitivity / (1 - specificity)
lr_negative LR- = (1 - sensitivity) / specificity
diagnostic_odds_ratio DOR = LR+ / LR-

Applies when

0 < specificity < 1
sensitivity < 1

Example

m.ask("Find the likelihood ratios of a test with sensitivity 0.9 and specificity 0.8.")

→ lr_positive = 4.5, lr_negative = 0.125, diagnostic_odds_ratio = 36

Odds and log-odds of a probabilityprobability_to_odds

The odds p/(1-p), the odds against and the log-odds (logit) of a probability.

\text{odds} = \frac{p}{1-p},\quad \operatorname{logit} p = \ln\frac{p}{1-p}

Inputs

probability probability p, strictly between 0 and 1

Returns

odds odds in favour p/(1-p)
odds_against odds against (1-p)/p
log_odds natural log-odds ln(p/(1-p))

Applies when

0 < p < 1

Example

m.ask("Convert the probability 0.2 to odds.")

→ odds = 0.25, odds_against = 4, log_odds = -1.38629

Post-test probability from a likelihood ratioposterior_from_likelihood_ratio

Bayes' theorem in odds form: post-test odds = pre-test odds x likelihood ratio, and the post-test probability.

\frac{P(H\mid E)}{1-P(H\mid E)} = \frac{P(H)}{1-P(H)}\times LR

Inputs

prior pre-test (prior) probability, 0 < p < 1
likelihood_ratio likelihood ratio (>= 0)

Returns

posterior post-test (posterior) probability
prior_odds pre-test odds
posterior_odds post-test odds

Applies when

0 < prior < 1
likelihood ratio >= 0

Example

m.ask("The pre-test probability is 0.2 and the likelihood ratio is 6. What is the post-test probability?")

→ posterior = 0.6, prior_odds = 0.25, posterior_odds = 1.5

Probability from oddsodds_to_probability

The probability o/(1+o) that corresponds to odds o in favour.

p = \frac{o}{1+o}

Inputs

odds odds in favour (>= 0)

Returns

probability probability o/(1+o)
complement 1 - probability

Applies when

odds >= 0

Example

m.ask("Convert odds of 3 to a probability.")

→ probability = 0.75, complement = 0.25

Posterior approximation 1

Laplace (normal) approximation of a beta posteriorlaplace_approximation_beta

Approximate a Beta(alpha, beta) posterior by a normal distribution centred at its mode with variance -1/(second derivative of the log density); the approximate credible interval.

\hat p = \frac{\alpha-1}{\alpha+\beta-2},\quad \sigma^2 = \left[\frac{\alpha-1}{\hat p^2} + \frac{\beta-1}{(1-\hat p)^2}\right]^{-1}

Inputs

alpha Beta alpha (> 1)
beta Beta beta (> 1)
level credibility level as a fraction (0.95 for 95%; default 0.95) optional

Returns

mode posterior mode
laplace_sd Laplace standard deviation
lower approximate lower limit
upper approximate upper limit

Applies when

alpha > 1 and beta > 1

Example

m.ask("Laplace approximation of a Beta(3, 3) posterior with a 95% interval.")

→ mode = 0.5, laplace_sd = 0.25, lower = 0.010009, upper = 0.989991

Prior elicitation 1

Beta prior from a mean and standard deviationbeta_prior_from_mean_sd

The Beta(alpha, beta) prior with a stated mean and standard deviation (method of moments) and its effective prior sample size alpha + beta.

\nu = \frac{m(1-m)}{s^2} - 1,\quad \alpha = m\nu,\quad \beta = (1-m)\nu

Inputs

mean prior mean, 0 < m < 1
sd prior standard deviation, sd^2 < m(1 - m)

Returns

alpha alpha
beta beta
prior_sample_size alpha + beta

Applies when

0 < mean < 1
0 < sd^2 < mean (1 - mean)

Example

m.ask("Find the beta prior with mean = 0.3 and standard deviation = 0.1.")

→ alpha = 6, beta = 14, prior_sample_size = 20

Sensitivity analysis 5

Morris elementary-effects screeningmorris_elementary_effects

Mean, mean absolute value (mu*) and standard deviation of the elementary effects from one-at-a-time trajectories: in each row, output j+1 is the model at the previous point with input j raised by delta.

EE_j = \frac{f(x + \Delta e_j) - f(x)}{\Delta},\quad \mu_j^* = \frac1r\sum_t |EE_j^{(t)}|

Inputs

outputs model outputs: one row per trajectory, k + 1 values (base point, then input 1, 2, ... raised in turn) matrix
delta step delta (> 0)

Returns

mu mean elementary effect of each input
mu_star mean absolute elementary effect
sigma standard deviation of the elementary effects

Applies when

at least 2 trajectories
delta > 0

Example

m.ask("Elementary effects screening with model outputs [[10, 12], [11, 10]] and step 1.")

→ mu = [0.5], mu_star = [1.5], sigma = [2.12132034356]

Sobol indices of a linear (additive) modelsobol_indices_linear_model

First-order Sobol indices of y = sum a_i x_i with independent inputs of given standard deviations (S_i = a_i^2 sigma_i^2 / Var y; they sum to 1).

S_i = \frac{a_i^2\sigma_i^2}{\sum_j a_j^2\sigma_j^2}

Inputs

coefficients model coefficients a_i vector
sds input standard deviations sigma_i (> 0) vector

Returns

first_order first-order (= total) indices
variance Var y

Applies when

independent inputs
sigma_i > 0

Example

m.ask("Sobol indices of a linear model with coefficients = [2, 1, 0.5] and sds = [1, 2, 4].")

→ first_order = [0.333333333333, 0.333333333333, 0.333333333333], variance = 12

Sobol indices of a model with an interaction termsobol_indices_interaction_model

First-order, interaction and total Sobol indices of y = a1 x1 + a2 x2 + a12 x1 x2 with independent normal inputs x1 ~ N(mean1, sd1^2), x2 ~ N(mean2, sd2^2).

V_1 = (a_1 + a_{12}\mu_2)^2\sigma_1^2,\ V_2 = (a_2 + a_{12}\mu_1)^2\sigma_2^2,\ V_{12} = a_{12}^2\sigma_1^2\sigma_2^2

Inputs

a1 coefficient of x1
a2 coefficient of x2
a12 interaction coefficient
mean1 mean of x1
sd1 sd of x1 (> 0)
mean2 mean of x2
sd2 sd of x2 (> 0)

Returns

S1 first-order index of x1
S2 first-order index of x2
S12 interaction index
ST1 total index of x1
ST2 total index of x2
variance Var y

Applies when

independent normal inputs

Example

m.ask("Sobol indices of the interaction model with a1 = 0, a2 = 0, a12 = 1, mean1 = 0, sd1 = 1, mean2 = 0, sd2 = 1.")

→ S1 = 0, S2 = 0, S12 = 1, ST1 = 1, ST2 = 1, variance = 1

Sobol indices of Sobol's G-functionsobol_indices_g_function

Exact first-order and total Sobol indices and the variance of the G-function prod_i (|4 x_i - 2| + a_i)/(1 + a_i), x_i uniform on [0, 1].

V_i = \frac{1}{3(1+a_i)^2},\ V = \prod_i(1+V_i) - 1,\ S_i = \frac{V_i}{V},\ S_{Ti} = \frac{V_i\prod_{j\ne i}(1+V_j)}{V}

Inputs

a parameters a_i >= 0 (small a_i = important input) vector

Returns

first_order first-order indices S_i
total_order total indices ST_i
variance variance of the function

Applies when

a_i >= 0
independent uniform inputs

Example

m.ask("Sobol indices of the G-function with a = [0, 0].")

→ first_order = [0.428571428571, 0.428571428571], total_order = [0.571428571429, 0.571428571429], variance = 0.777778

Sobol sensitivity indices of the Ishigami functionsobol_indices_ishigami

Exact first-order and total Sobol indices and the variance of the Ishigami test function f = sin x1 + a sin^2 x2 + b x3^4 sin x1 with x_i uniform on [-pi, pi].

V = \frac{a^2}{8} + \frac{b\pi^4}{5} + \frac{b^2\pi^8}{18} + \frac12,\ V_1 = \frac12\left(1 + \frac{b\pi^4}{5}\right)^2,\ V_2 = \frac{a^2}{8},\ V_{13} = \frac{8b^2\pi^8}{225}

Inputs

a parameter a (commonly 7)
b parameter b (commonly 0.1)

Returns

S1 first-order index of x1
S2 first-order index of x2
S3 first-order index of x3
ST1 total index of x1
ST2 total index of x2
ST3 total index of x3
variance variance of f

Applies when

inputs independent uniform on [-pi, pi]

Example

m.ask("Sobol indices of the Ishigami function with a = 7 and b = 0.1.")

→ S1 = 0.313905, S2 = 0.442411, S3 = 0, ST1 = 0.557589, ST2 = 0.442411, ST3 = 0.243684, variance = 13.8446

Uncertainty propagation 14

Covariance propagation through a linear map y = A xcovariance_propagation_linear

Covariance matrix of the outputs y = A x and their standard uncertainties, given the covariance matrix of x (first-order propagation with a sensitivity matrix).

\Sigma_y = A\,\Sigma_x\,A^{\mathsf T}

Inputs

matrix sensitivity (Jacobian) matrix A: one row per output matrix
covariance covariance matrix of the inputs (symmetric, positive semi-definite) matrix

Returns

output_covariance A Sigma A^T
uncertainties square roots of its diagonal

Applies when

covariance symmetric positive semi-definite
columns of A = inputs

Example

m.ask("Covariance propagation: matrix = [[1, 1], [1, -1]], covariance = [[0.04, 0.01], [0.01, 0.09]].")

→ output_covariance = [[0.15, -0.05], [-0.05, 0.11]], uncertainties = [0.387298334621, 0.331662479036]

Effective degrees of freedom and expanded uncertainty (Welch-Satterthwaite)welch_satterthwaite_expanded_uncertainty

Combined standard uncertainty, Welch-Satterthwaite effective degrees of freedom, the Student-t coverage factor and the expanded uncertainty at a stated coverage probability.

\nu_{eff} = \frac{u_c^4}{\sum_i u_i^4/\nu_i},\quad U = t_{(1+p)/2,\,\nu_{eff}}\,u_c

Inputs

uncertainties uncertainty components c_i u_i (> 0) vector
dof degrees of freedom of each component (> 0) vector
level credibility level as a fraction (0.95 for 95%; default 0.95) optional

Returns

combined_uncertainty u_c
effective_dof nu_eff
coverage_factor k = t_{(1+p)/2, nu_eff}
expanded_uncertainty U = k u_c

Applies when

components and degrees of freedom > 0
0 < level < 1

Example

m.ask("Welch-Satterthwaite: uncertainties = [1], degrees of freedom = [10]. Expanded uncertainty at 95% coverage?")

→ combined_uncertainty = 1, effective_dof = 10, coverage_factor = 2.22814, expanded_uncertainty = 2.22814

Exact mean and standard deviation of a product of independent normalsproduct_of_normals_exact_moments

Exact moments of XY for independent normal X and Y, compared with first-order propagation (which misses the sd_x sd_y term).

\mathrm{Var}(XY) = \mu_x^2\sigma_y^2 + \mu_y^2\sigma_x^2 + \sigma_x^2\sigma_y^2

Inputs

mean_x mean of X
sd_x sd of X (> 0)
mean_y mean of Y
sd_y sd of Y (> 0)

Returns

exact_mean mean_x mean_y
exact_sd sqrt(mx^2 sy^2 + my^2 sx^2 + sx^2 sy^2)
first_order_sd sqrt(mx^2 sy^2 + my^2 sx^2)

Applies when

independent normal inputs

Example

m.ask("Exact moments of the product of independent normals: mean_x = 10, sd_x = 1, mean_y = 5, sd_y = 2.")

→ exact_mean = 50, exact_sd = 20.7123, first_order_sd = 20.6155

Exact mean and variance of a polynomial of a normal variablepolynomial_of_normal_moments

Exact mean, variance and standard deviation of Y = c0 + c1 X + ... + c6 X^6 for X ~ N(mean, sd^2) (nonlinear uncertainty propagation without linearization).

E[Z^{2m}] = (2m-1)!!,\quad X = \mu + \sigma Z

Inputs

coefficients coefficients c0, c1, ... (constant first, up to degree 6) vector
mean mean of X
sd standard deviation of X (> 0)

Returns

mean E[Y]
variance Var[Y]
sd sd of Y

Applies when

sd > 0
degree <= 6

Example

m.ask("Mean and variance of the polynomial with coefficients = [1, 2, 3] of a normal X with mean = 0, sd = 1.")

→ mean = 4, variance = 22, sd = 4.69042

Exact versus first-order propagation through exp(X)lognormal_propagation

Mean and standard deviation of Y = exp(X) for normal X ~ N(mu, sigma^2), exactly (lognormal) and by first-order propagation, showing the bias of the linear approximation.

E[e^X] = e^{\mu + \sigma^2/2},\quad \mathrm{SD}[e^X] = e^{\mu + \sigma^2/2}\sqrt{e^{\sigma^2} - 1}

Inputs

mu mean of X
sigma standard deviation of X (> 0)

Returns

exact_mean exp(mu + sigma^2/2)
exact_sd sqrt(exp(sigma^2) - 1) exp(mu + sigma^2/2)
first_order_mean exp(mu)
first_order_sd exp(mu) sigma

Applies when

sigma > 0

Example

m.ask("Lognormal propagation with mu = 0 and sigma = 1.")

→ exact_mean = 1.64872, exact_sd = 2.1612, first_order_mean = 1, first_order_sd = 1

Inverse-variance weighted mean of measurementsinverse_variance_weighted_mean

Combine measurements of one quantity with different standard uncertainties: weighted mean, its uncertainty, chi-square of consistency and Birge ratio (equal to the posterior under a flat prior).

\bar x_w = \frac{\sum x_i/u_i^2}{\sum 1/u_i^2},\quad u(\bar x_w) = \Big(\sum 1/u_i^2\Big)^{-1/2}

Inputs

values measured values vector
uncertainties their standard uncertainties (> 0) vector

Returns

weighted_mean sum(x_i/u_i^2)/sum(1/u_i^2)
uncertainty 1/sqrt(sum 1/u_i^2)
chi_square sum((x_i - mean)^2/u_i^2)
birge_ratio sqrt(chi^2/(n-1))

Applies when

u_i > 0
n >= 2

Example

m.ask("Combine the measurements [5.1, 4.9] with uncertainties [0.1, 0.1]: weighted mean?")

→ weighted_mean = 5, uncertainty = 0.0707107, chi_square = 2, birge_ratio = 1.41421

Type B standard uncertainty from a stated half-widthtype_b_standard_uncertainty

Standard uncertainty of a quantity known only to lie within +/- a, for a rectangular, triangular or U-shaped (arcsine) distribution.

u = \frac{a}{\sqrt 3}\ (\text{rect.}),\ \frac{a}{\sqrt 6}\ (\text{tri.}),\ \frac{a}{\sqrt 2}\ (\text{arcsine})

Inputs

half_width half-width a of the interval (> 0)
distribution rectangular, triangular or arcsine (U-shaped) choice

Returns

standard_uncertainty a / divisor
divisor sqrt(3), sqrt(6) or sqrt(2)

Applies when

a > 0

Example

m.ask("Type B uncertainty: half-width = 0.5 with a rectangular distribution.")

→ standard_uncertainty = 0.288675, divisor = 1.73205

Uncertainty of a function of one variableuncertainty_function

Value and first-order standard uncertainty of f(x) for exp, ln, log10, sqrt, sin, cos, tan (x in radians) or 1/x.

u_f = \left|f'(x)\right| u_x

Inputs

x value of x (radians for sin, cos, tan)
ux standard uncertainty of x
function exp, ln, log10, sqrt, sin, cos, tan or reciprocal choice

Returns

value f(x)
uncertainty |f'(x)| u(x)

Applies when

x inside the function's domain

Example

m.ask("Uncertainty of a function: function = ln, x = 2.5, ux = 0.05.")

→ value = 0.916291, uncertainty = 0.02

Uncertainty of a linear combination of correlated inputsuncertainty_linear_correlated

Combined standard uncertainty of y = sum c_i x_i when the inputs are correlated (correlation matrix), u_y^2 = sum_ij c_i c_j u_i u_j r_ij.

u_y^2 = \sum_i\sum_j c_i c_j u_i u_j r_{ij}

Inputs

values input values vector
uncertainties standard uncertainties vector
correlation correlation matrix (ones on the diagonal, positive semi-definite) matrix
coefficients coefficients c_i (default 1) vector optional

Returns

value y
uncertainty u_y

Applies when

valid correlation matrix

Example

m.ask("Uncertainty of a linear combination with correlation matrix [[1, 0.5], [0.5, 1]], values = [10, 20], uncertainties = [0.3, 0.4], coefficients = [1, 1].")

→ value = 30, uncertainty = 0.608276

Uncertainty of a power x^nuncertainty_power

Value, standard uncertainty and relative uncertainty of y = x^n (first order).

\frac{u_y}{|y|} = |n|\,\frac{u_x}{|x|}

Inputs

x value of x
ux standard uncertainty of x (>= 0)
exponent exponent n

Returns

value x^n
uncertainty |n x^(n-1)| u(x)
relative_uncertainty |n| u(x)/|x|

Applies when

x > 0 for a non-integer exponent
x != 0 for n < 1

Example

m.ask("Uncertainty of a power: x = 3.0, ux = 0.1, exponent = 2.")

→ value = 9, uncertainty = 0.6, relative_uncertainty = 0.0666667

Uncertainty of a power-law product y = k x1^p1 x2^p2 ...uncertainty_power_law

Value and first-order uncertainty of a product of powers of independent positive inputs (relative uncertainties weighted by the exponents add in quadrature).

\frac{u_y}{|y|} = \sqrt{\sum_i \left(p_i\,\frac{u_i}{x_i}\right)^2}

Inputs

coefficient exact constant k
values positive input values x_i vector
exponents exponents p_i vector
uncertainties standard uncertainties u_i vector

Returns

value y
uncertainty u_y
relative_uncertainty u_y/|y|

Applies when

x_i > 0
independent inputs

Example

m.ask("Uncertainty of a power law: coefficient = 0.5, values = [2.0, 3.0], exponents = [1, 2], uncertainties = [0.02, 0.06].")

→ value = 9, uncertainty = 0.37108, relative_uncertainty = 0.0412311

Uncertainty of a product and a quotientuncertainty_product_quotient

Value and first-order standard uncertainty of a*b and a/b for independent a and b.

\frac{u_{ab}}{|ab|} = \frac{u_{a/b}}{|a/b|} = \sqrt{\left(\frac{u_a}{a}\right)^2 + \left(\frac{u_b}{b}\right)^2}

Inputs

a value of a
ua standard uncertainty of a (>= 0)
b value of b (non-zero)
ub standard uncertainty of b (>= 0)

Returns

product a*b
u_product u(a*b)
quotient a/b
u_quotient u(a/b)

Applies when

independent a and b
b != 0

Example

m.ask("Uncertainty of the product and quotient: a = 12.0, ua = 0.3, b = 4.0, ub = 0.2.")

→ product = 48, u_product = 2.68328, quotient = 3, u_quotient = 0.167705

Uncertainty of a product and quotient with correlationuncertainty_product_correlated

First-order standard uncertainty of a*b and a/b when a and b are correlated with coefficient rho.

u_{ab}^2 = b^2u_a^2 + a^2u_b^2 + 2ab\,\rho\,u_au_b

Inputs

a value of a
ua standard uncertainty of a
b value of b (non-zero)
ub standard uncertainty of b
rho correlation coefficient of a and b (-1 to 1)

Returns

product a*b
u_product u(a*b)
quotient a/b
u_quotient u(a/b)

Applies when

-1 <= rho <= 1
b != 0

Example

m.ask("Uncertainty of the product with correlation: a = 1, ua = 0.1, b = 1, ub = 0.1, rho = -1.")

→ product = 1, u_product = 0, quotient = 1, u_quotient = 0.2

Uncertainty of a sum or linear combination (independent inputs)uncertainty_sum

Value and combined standard uncertainty of y = sum c_i x_i for independent inputs (default c_i = 1; use -1 for a difference).

u_y = \sqrt{\sum_i c_i^2 u_i^2}

Inputs

values input values x_i vector
uncertainties standard uncertainties u_i (>= 0) vector
coefficients coefficients c_i (default all 1) vector optional

Returns

value y = sum c_i x_i
uncertainty u_y = sqrt(sum (c_i u_i)^2)

Applies when

independent inputs
u_i >= 0

Example

m.ask("Uncertainty of a sum: values = [12.3, 4.5, 7.1], uncertainties = [0.2, 0.1, 0.3].")

→ value = 23.9, uncertainty = 0.374166

Versions & files

Releases and repository files

Files

15 files, 36.2 MB in total. The largest:

File	Size
model.safetensors	18.9 MB
dist/bayesian_statistics_ai999x-1.0.0-py3-none-any.whl	8.7 MB
onnx/search/encoder.int8.onnx	4.9 MB
docs/bayesian_statistics_ai999x_capabilities.pdf	2.7 MB
onnx/search/index.json	403.3 KB
docs/bayesian_statistics_ai999x_capabilities.html	236.2 KB
docs/bayesian_statistics_ai999x_capabilities.md	178.7 KB
onnx/search/tokenizer.json	176.4 KB
README.md	22.8 KB
LICENSE	11.4 KB
evaluation/release_gate.md	2.3 KB
evaluation/KNOWLEDGE_CARD.md	1.4 KB

Browse every file ↗

Questions

Frequently asked

What does bayesian-statistics-ai999x do?

Bayesian Statistics answers questions in bayesian statistics: Bayes' theorem, conjugate priors and posteriors, credible intervals, Bayes factors, Monte Carlo estimation, uncertainty propagation, sensitivity analysis and decision analysis. Each answer gives the value with its unit, the steps and the formula in LaTeX.

How do I install bayesian-statistics-ai999x?

With pip: pip install bayesian-statistics-ai999x. The same files are on Hugging Face (hf download AI999X/bayesian-statistics-ai999x --local-dir bayesian-statistics-ai999x) and GitHub.

Does it run offline?

Yes. It is a Python package that runs on an ordinary CPU without a network connection (Python >=3.10); no GPU and no API key are needed.

How many capabilities does bayesian-statistics-ai999x have?

78 capabilities in 12 topics, including Conjugate priors, Uncertainty propagation, Decision analysis, Credible intervals, Bayes factors and model comparison.

Which field and tasks does it cover?

Field: Mathematics & Statistics. Tasks: Question Answering, Sentence Similarity, Feature Extraction.

What does an answer look like?

Asked "A test has sensitivity 0.9 and specificity 0.95; the prevalence is 0.01. What is the positive predictive value?", it answers: ppv = 0.153846, npv = 0.998938, p_positive = 0.0585, false_discovery_rate = 0.846154, false_omission_rate = 0.00106213.

What happens when a question cannot be answered?

It says so and names what is missing or out of range, instead of guessing a number.

Explore every model.

More in Mathematics & Statistics All models Capability orb