physics.statistics

physics.statistics

Statistical and workflow propagation helpers.

Functions

Name	Description
calc_arcsine_two_arm_sample_size	Plan equal arm sizes using an arcsine effect size.
calc_behavior_cloning_cost_bound	Evaluate quadratic horizon growth of imitation excess cost.
calc_cdf_quantile	Select the first value attaining a supplied cumulative fraction.
calc_chance_agreement	Compute chance agreement from a shared label distribution.
calc_chi2_gate	Compute a chi-square consistency-test threshold.
calc_chi2_mean_band	Bound the mean of independent chi-square consistency statistics.
calc_chi_square_quantile	Compute a Gaussian confidence region’s squared Mahalanobis radius.
calc_clopper_pearson_interval	Compute a two-sided exact binomial confidence interval.
calc_cohens_h	Compute the signed arcsine effect size for two proportions.
calc_cohens_kappa	Normalize observed agreement relative to chance agreement.
calc_conjunctive_confidence	Combine independent premise probabilities by conjunction.
calc_consecutive_detection_frames	Find the frame count for an all-miss probability target.
calc_cross_entropy_exponential_proposal	Trace population cross-entropy updates for an exponential tail.
calc_dagger_cost_bound	Evaluate linear horizon growth under recoverable imitation.
calc_effective_sample_size	Approximate the independent sample count of an AR(1) mean.
calc_equal_n_for_z	Solve the equal-group size for a specified pooled z statistic.
calc_expected_calibration_error	Compute count-weighted absolute calibration error.
calc_exponential_decay_rate	Fit an exponential rate through two readings.
calc_fisher_exact_two_sided	Compute a two-sided Fisher exact contingency-table p-value.
calc_gaussian_kl	Compute the KL divergence between scalar normal distributions.
calc_gaussian_kl_terms	Decompose the KL divergence between scalar normal distributions.
calc_importance_sampling_exponential	Compute exponential-tail importance-sampling moments.
calc_jensen_shannon_divergence	Measure divergence between aligned probability distributions.
calc_ks_coefficient	Compute an asymptotic two-sample KS critical coefficient.
calc_ks_critical_value	Scale a two-sample Kolmogorov-Smirnov critical coefficient.
calc_ks_normal_location_shift	Compute the CDF distance for an equal-variance normal shift.
calc_linf_corner_l2	Convert an infinity-norm corner radius to its Euclidean norm.
calc_mcnemar_exact_two_sided	Compute the exact conditional McNemar p-value.
calc_mcnemar_max_discordant	Bound discordance for an uncorrected McNemar threshold.
calc_mean_half_width	Calculate a Student t interval half-width for a mean.
calc_monte_carlo_sample_size	Estimate independent draws for a relative standard-error target.
calc_multilegged_confidence	Combine independent supporting legs by disjunction.
calc_newcombe_paired_difference	Estimate a paired rate difference and its interval.
calc_noncentral_gate_detection	Compute the probability that a biased statistic crosses a gate.
calc_one_proportion_sample_size	Plan a one-sample binomial comparison of two rates.
calc_optimal_exponential_proposal_mean	Minimize the exponential importance-sampling second moment.
calc_paired_difference_variance	Calculate variance of a signed paired outcome.
calc_poisson_exact_interval	Invert Poisson tails to bound an event-count mean.
calc_poisson_mean_for_reserve	Find a mean-count boundary for a fixed Poisson reserve.
calc_poisson_reserve	Find the smallest reserve satisfying a Poisson tail target.
calc_poisson_tail	Compute the probability of at least a specified event count.
calc_population_stability_index	Measure drift between two aligned binned distributions.
calc_proportion_gap_interval	Bound a proportion difference using its unpooled standard error.
calc_proportion_gap_se	Estimate the standard error of two independent proportions.
calc_proportion_half_width	Estimate the normal interval half-width for a proportion.
calc_sample_collection_time	Estimate collection time for a fixed traffic allocation.
calc_samples_for_half_width	Plan independent Bernoulli draws for a target half-width.
calc_samples_per_arm_for_difference	Plan equal-arm observations for a mean difference.
calc_smoothing_radius	Compute the binary Gaussian randomized-smoothing radius.
calc_spearman_rank_correlation	Correlate paired observations through their average ranks.
calc_split_conformal_min_count	Find the minimum calibration count for a finite conformal rank.
calc_split_conformal_rank	Find the finite-sample split-conformal calibration rank.
calc_sprt_boundaries	Compute Wald stopping boundaries for a log-likelihood ratio.
calc_sprt_expected_samples	Approximate a likelihood walk stopping count.
calc_sprt_increments	Compute Bernoulli SPRT log-likelihood increments.
calc_standard_error	Compute the standard error of an independent sample mean.
calc_t_critical	Return a two-sided Student t critical value.
calc_task_averaged_variance	Calculate variance of an equally weighted task average.
calc_two_proportion_detectable_lift	Invert the common-variance two-arm planning formula.
calc_two_proportion_sample_size	Estimate per-arm sample size for a two-proportion comparison.
calc_two_proportion_z	Compute a pooled two-proportion z statistic.
calc_wald_interval	Compute an unclipped Wald interval for a binomial proportion.
calc_welch_interval	Construct a Welch interval for two independent means.
calc_wilson_clean_runs	Find the clean-run count for a Wilson upper bound.
calc_wilson_interval	Compute a Wilson score interval for a binomial proportion.
calc_within_task_variance	Recover within-task Bernoulli variance from task heterogeneity.
fit_exponential_growth	Fit a constant multiplicative growth trend in log space.

calc_arcsine_two_arm_sample_size

physics.statistics.calc_arcsine_two_arm_sample_size(
    effect_h,
    alpha=0.05,
    power=0.8,
)

Plan equal arm sizes using an arcsine effect size.

Parameters

Name	Type	Description	Default
effect_h	float	Nonzero dimensionless Cohen h effect size; its sign does not affect sample size.	required
alpha	float	Two-sided significance probability in (0, 1); default 0.05.	0.05
power	float	Target power probability in (0, 1); default 0.80.	0.8

Returns

Name	Type	Description
	int	Required observations per arm, rounded upward.

Notes

Equation: n = ceil(2*((Phi**(-1)(1-alpha/2)+Phi**(-1)(power))/h)**2).

Assumptions: Independent equal-size arms and the large-sample variance approximation for transformed proportions.

Excludes: Exact binomial power, unequal allocation, multiple testing, and input-domain validation.

References

Cohen (1988), Statistical Power Analysis for the Behavioral Sciences, chapter 6. Derivation: Each transformed rate has variance 1/n; solve the two-arm normal power relation with variance 2/n.

Examples

>>> calc_arcsine_two_arm_sample_size(0.4) < calc_arcsine_two_arm_sample_size(0.2)
True

calc_behavior_cloning_cost_bound

physics.statistics.calc_behavior_cloning_cost_bound(
    horizon_steps,
    per_step_error,
)

Evaluate quadratic horizon growth of imitation excess cost.

Parameters

Name	Type	Description	Default
horizon_steps	int or float	Nonnegative finite-horizon decision count S.	required
per_step_error	float	Demonstrator-distribution classification error epsilon in [0, 1].	required

Returns

Name	Type	Description
	float	Worst-case excess cumulative cost under normalized per-step costs.

Notes

Equation: Excess cost bound = S**2*epsilon.

Assumptions: Supervised imitation on expert states, bounded per-step cost, and possible persistent state-distribution shift after an error.

Excludes: Task-specific recovery, learned-policy guarantees, and finite-data estimation of epsilon.

Ross and Bagnell (2010), Theorem 2.1: a policy trained only on the demonstrator’s states compounds a per-step error eps quadratically.

References

Ross and Bagnell (2010), Efficient Reductions for Imitation Learning. https://proceedings.mlr.press/v9/ross10a.html

Examples

>>> calc_behavior_cloning_cost_bound(20, 0.01) == 4 * calc_behavior_cloning_cost_bound(10, 0.01)
True

calc_cdf_quantile

physics.statistics.calc_cdf_quantile(values, cumulative, fraction)

Select the first value attaining a supplied cumulative fraction.

Parameters

Name	Type	Description	Default
values	sequence	Nonempty support values in ascending order and common units.	required
cumulative	sequence of float	Aligned nondecreasing cumulative probabilities in [0, 1].	required
fraction	float	Requested quantile fraction in [0, 1].	required

Returns

Name	Type	Description
	object	First support value at or above the fraction; last value when the supplied CDF never reaches it.

Notes

Equation: Return values[min{i : cumulative[i]>=fraction}], or the final value if no crossing exists.

Assumptions: Caller supplies aligned ordered CDF data; no interpolation is performed.

Excludes: CDF estimation, input normalization, interpolation, and uncertainty intervals.

values ascend with cumulative (a released CDF); past the end the last value is returned.

References

Derivation: the generalized inverse of a discrete CDF selects its first threshold crossing; endpoint fallback is the implementation convention.

Examples

>>> calc_cdf_quantile([1, 3], [0.5, 1.0], 0.75)
3

calc_chance_agreement

physics.statistics.calc_chance_agreement(class_shares)

Compute chance agreement from a shared label distribution.

Parameters

Name	Type	Description	Default
class_shares	sequence of float	Nonempty nonnegative class probabilities summing to one.	required

Returns

Name	Type	Description
	float	Dimensionless independent chance-agreement probability.

Notes

Equation: p_e = sum(s_k**2).

Assumptions: Annotators independently sample labels from the same marginal distribution.

Excludes: Distinct annotator marginals, label dependence, and estimation uncertainty.

class_shares is each annotator’s label mix (fractions summing to 1); both annotators are assumed to share it.

References

Derivation: matching on class k has probability s_k*s_k; sum mutually exclusive matches.

Examples

>>> calc_chance_agreement([0.5, 0.5])
0.5

calc_chi2_gate

physics.statistics.calc_chi2_gate(confidence, dof)

Compute a chi-square consistency-test threshold.

Parameters

Name	Type	Description	Default
confidence	float	Desired lower-tail probability in (0, 1).	required
dof	int	Positive number of degrees of freedom.	required

Returns

Name	Type	Description
	float	Dimensionless chi-square threshold; values above it fail the gate.

Notes

Equation: gate = F_chi2(dof)**(-1)(confidence).

Assumptions: The normalized squared statistic has a central chi-square distribution with the specified degrees of freedom.

Excludes: Verification of Gaussian residuals or covariance calibration, repeated-testing corrections, and bias modeling.

References

Bar-Shalom, Li, and Kirubarajan (2001), Estimation with Applications to Tracking and Navigation, section 5.4. SciPy chi-square distribution and inverse CDF: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.chi2.html

Examples

>>> calc_chi2_gate(0.99, 2) > calc_chi2_gate(0.95, 2)
True

calc_chi2_mean_band

physics.statistics.calc_chi2_mean_band(n_samples, dof, confidence)

Bound the mean of independent chi-square consistency statistics.

Parameters

Name	Type	Description	Default
n_samples	int	Positive number of independent consistency statistics.	required
dof	int	Positive degrees of freedom of each statistic.	required
confidence	float	Central interval probability in (0, 1).	required

Returns

Name	Type	Description
	tuple of float	Dimensionless lower and upper thresholds for the sample mean.

Notes

Equation: band = (chi2_ppf((1-c)/2, N*d), chi2_ppf((1+c)/2, N*d))/N.

Assumptions: Independent central chi-square statistics with a common dimension d; their sum has N*d degrees of freedom.

Excludes: Serially correlated residuals, estimated effective sample size, and noncentral alternatives.

References

Bar-Shalom, Li, and Kirubarajan (2001), section 5.4.2. SciPy chi-square distribution and inverse CDF: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.chi2.html

Examples

>>> low, high = calc_chi2_mean_band(100, 2, 0.95)
>>> low < 2 < high
True

calc_chi_square_quantile

physics.statistics.calc_chi_square_quantile(degrees_of_freedom, coverage)

Compute a Gaussian confidence region’s squared Mahalanobis radius.

Parameters

Name	Type	Description	Default
degrees_of_freedom	int	Positive Gaussian dimension.	required
coverage	float	Probability enclosed by the confidence region, in (0, 1).	required

Returns

Name	Type	Description
	float	Dimensionless squared Mahalanobis radius; the physical semiaxes are sqrt(radius_squared) times covariance-axis standard deviations.

Notes

Equation: r_squared = F_chi2(k)**(-1)(coverage).

Assumptions: A Gaussian vector with known mean and nonsingular covariance; the whitened squared norm is chi-square.

Excludes: Uncertainty from covariance estimation, non-Gaussian coverage, and physical axis conversion.

References

SciPy chi-square distribution and inverse CDF: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.chi2.html

Examples

>>> calc_chi_square_quantile(2, 0.95) == calc_chi2_gate(0.95, 2)
True

calc_clopper_pearson_interval

physics.statistics.calc_clopper_pearson_interval(
    successes,
    trials,
    confidence=0.95,
)

Compute a two-sided exact binomial confidence interval.

Exact two-sided Clopper-Pearson interval for successes out of trials.

Equation (beta quantiles, alpha = 1 - confidence): lower = B^{-1}(alpha/2; k, n - k + 1) (0 when k = 0) upper = B^{-1}(1 - alpha/2; k + 1, n - k) (1 when k = n)

Source: Clopper and Pearson, “The Use of Confidence or Fiducial Limits Illustrated in the Case of the Binomial,” Biometrika 26(4), 1934.

Parameters

Name	Type	Description	Default
successes	float	Integer successes from zero through trials.	required
trials	float	Positive integer trial count.	required
confidence	float	Central confidence probability in (0, 1).	0.95

Returns

Name	Type	Description
	tuple of float	Lower and upper binomial probability limits.

Notes

Equation: Limits are beta quantiles at alpha/2 and 1-alpha/2 with endpoint conventions.

Assumptions: Independent identically distributed Bernoulli trials and fixed sample size.

Excludes: Sequential testing and dependent observations.

References

Derivation: Clopper and Pearson (1934), The Use of Confidence or Fiducial Limits Illustrated in the Case of the Binomial.

Examples

>>> calc_clopper_pearson_interval(0, 10)[0] == 0
True

calc_cohens_h

physics.statistics.calc_cohens_h(p_high, p_low)

Compute the signed arcsine effect size for two proportions.

Parameters

Name	Type	Description	Default
p_high	float	First probability on [0, 1]; the name does not enforce ordering.	required
p_low	float	Second probability on [0, 1].	required

Returns

Name	Type	Description
	float	Dimensionless signed Cohen h; positive when p_high exceeds p_low.

Notes

Equation: h = 2*asin(sqrt(p_high)) - 2*asin(sqrt(p_low)).

Assumptions: Probabilities describe comparable binary outcomes.

Excludes: Sampling uncertainty, a test decision, and validation of probability bounds.

References

Cohen (1988), Statistical Power Analysis for the Behavioral Sciences, chapter 6. Derivation: Subtract the two variance-stabilizing arcsine-transformed probabilities.

Examples

>>> calc_cohens_h(0.5, 0.5) == 0
True

calc_cohens_kappa

physics.statistics.calc_cohens_kappa(observed_agreement, chance_agreement)

Normalize observed agreement relative to chance agreement.

Parameters

Name	Type	Description	Default
observed_agreement	float	Observed matching-label fraction in [0, 1].	required
chance_agreement	float	Chance agreement probability in [0, 1).	required

Returns

Name	Type	Description
	float	Dimensionless kappa; negative values indicate agreement below chance.

Notes

Equation: kappa = (p_o-p_e)/(1-p_e).

Assumptions: Chance agreement is estimated from suitable annotator marginals.

Excludes: Confidence intervals, weighted ordinal disagreement, and prevalence adjustment.

References

Derivation: scale observed excess agreement p_o-p_e by the maximum possible excess 1-p_e. Cohen (1960) introduced this agreement coefficient.

Examples

>>> calc_cohens_kappa(0.75, 0.5)
0.5

calc_conjunctive_confidence

physics.statistics.calc_conjunctive_confidence(premise_confidences)

Combine independent premise probabilities by conjunction.

Parameters

Name	Type	Description	Default
premise_confidences	sequence of float	Nonempty probabilities p_i in [0, 1].	required

Returns

Name	Type	Description
	float	Dimensionless probability that every premise holds.

Notes

Equation: P(all premises) = product(p_i).

Assumptions: Independent premise events, each necessary for the conclusion.

Excludes: Common-cause dependence and validation of supplied confidence judgments.

Equation: P(claim) = prod_i p_i

The conjunction is no stronger than its weakest premise and, with more than one premise below 1, strictly weaker.

Source: probability of the intersection of independent events, as applied to conjunctive assurance arguments in Bloomfield and Littlewood 2003, “Multi-legged arguments: the impact of diversity upon confidence in dependability arguments” (DSN 2003).

References

Derivation: repeated multiplication of independent-event intersection probabilities.

Examples

>>> calc_conjunctive_confidence([0.5, 0.5])
0.25

calc_consecutive_detection_frames

physics.statistics.calc_consecutive_detection_frames(
    per_frame_detection,
    all_miss_probability,
)

Find the frame count for an all-miss probability target.

Consecutive frames a person must be in view before detection may be counted on.

Equation: N = ceil( ln p_miss / ln(1 - p_frame) )

the fewest frames for which the chance of missing the person in every one, (1 - p_frame)^N, is at most all_miss_probability, with misses independent from frame to frame (the best case; correlated misses need more frames). A tolerance of 1e-12 keeps an exact ratio from rounding up.

Parameters

Name	Type	Description	Default
per_frame_detection	float	Detection probability per independent frame, in (0, 1).	required
all_miss_probability	float	Maximum all-miss probability in (0, 1).	required

Returns

Name	Type	Description
	int	Minimum frame count under the independent-miss model.

Notes

Equation: N = ceil(log(all_miss_probability)/log(1-per_frame_detection)-1e-12).

Assumptions: Independent frame misses with a constant detection probability.

Excludes: Temporal correlation and uncertainty in detector performance.

References

Derivation: Multiply independent miss probabilities and invert the geometric decay.

Examples

>>> calc_consecutive_detection_frames(0.5, 0.25) == 2
True

calc_cross_entropy_exponential_proposal

physics.statistics.calc_cross_entropy_exponential_proposal(
    threshold,
    operational_mean,
    elite_share,
)

Trace population cross-entropy updates for an exponential tail.

Parameters

Name	Type	Description	Default
threshold	float or Quantity	Positive target failure threshold gamma.	required
operational_mean	float or Quantity	Positive operational mean m0, in threshold units.	required
elite_share	float	Elite proposal fraction beta in (0, 1).	required

Returns

Name	Type	Description
	dict	Levels contains thresholds in input units, n_levels is their count, and proposal_mean is the final mean in input units.

Notes

Equation: level_k = min(gamma, -m_k*log(beta)); m_(k+1) = level_k + m0, until level_k reaches gamma.

Assumptions: Population expectations replace empirical elites; the initial proposal mean equals the operational mean.

Excludes: Sampling variability, weight degeneracy, finite-sample stopping, and non-exponential families.

From q_0 = p = Exp(m0), each round sets the level to the (1 - beta) quantile of the current proposal, level = -m_k ln(beta) (capped at the threshold), and refits the proposal mean to the likelihood-weighted elite draws; for this family that is E_p[D | D > level] = level + m0. Stops when the level reaches the threshold, so the final mean is threshold + m0.

Returns a dict with levels (list, in the units of threshold), n_levels, and proposal_mean.

Source: Rubinstein and Kroese, The Cross-Entropy Method (Springer, 2004), multilevel rare-event algorithm; the exponential refit is in closed form.

References

Derivation: exponential memorylessness gives E[X | X > level] = level + m0; the proposal quantile supplies the next level. Rubinstein and Kroese (2004), The Cross-Entropy Method, motivates multilevel adaptation.

Examples

>>> result = calc_cross_entropy_exponential_proposal(2.0, 1.0, 0.5)
>>> result["levels"][-1] == 2.0 and result["proposal_mean"] == 3.0
True

calc_dagger_cost_bound

physics.statistics.calc_dagger_cost_bound(
    horizon_steps,
    per_step_error,
    recoverability=1.0,
)

Evaluate linear horizon growth under recoverable imitation.

Parameters

Name	Type	Description	Default
horizon_steps	int or float	Nonnegative decision horizon S.	required
per_step_error	float	Learner-state error epsilon in [0, 1].	required
recoverability	float	Nonnegative excess-to-go bound u per error; default 1.0.	1.0

Returns

Name	Type	Description
	float	Excess cumulative cost bound in the cost units of recoverability.

Notes

Equation: Excess cost bound = u*S*epsilon.

Assumptions: Expert labels on learner-visited states and a bounded recoverability constant u.

Excludes: Finite-iteration regret terms, policy optimization, and proof that a given system satisfies the assumptions.

Ross, Gordon, and Bagnell (2011), Theorem 3.2: linear in the horizon while the recoverability constant u stays bounded.

References

Ross, Gordon, and Bagnell (2011), A Reduction of Imitation Learning and Structured Prediction to No-Regret Online Learning. https://proceedings.mlr.press/v15/ross11a/ross11a.pdf

Examples

>>> calc_dagger_cost_bound(20, 0.01) == 2 * calc_dagger_cost_bound(10, 0.01)
True

calc_effective_sample_size

physics.statistics.calc_effective_sample_size(n_records, lag1_correlation)

Approximate the independent sample count of an AR(1) mean.

Parameters

Name	Type	Description	Default
n_records	int or float	Nonnegative number of observations.	required
lag1_correlation	float	Stationary lag-one correlation rho in (-1, 1).	required

Returns

Name	Type	Description
	float	Effective count under the asymptotic AR(1) variance approximation.

Notes

Equation: N_eff = N*(1-rho)/(1+rho).

Assumptions: A stationary AR(1) autocorrelation rho**k and enough records for the infinite-lag approximation.

Excludes: Finite-sample edge corrections, nonstationary drift, and arbitrary dependence.

N_eff = N (1 - rho) / (1 + rho) for estimating the mean of a series whose adjacent records have lag-one correlation rho (Thiebaux and Zwiers 1984).

References

Derivation: sum 1+2*sum(rho**k, k>=1)=(1+rho)/(1-rho) as the variance inflation of the sample mean; divide N by that inflation.

Examples

>>> calc_effective_sample_size(100, 0.0)
100.0

calc_equal_n_for_z

physics.statistics.calc_equal_n_for_z(rate_a, rate_b, z_critical)

Solve the equal-group size for a specified pooled z statistic.

Parameters

Name	Type	Description	Default
rate_a	float	Distinct probabilities on [0, 1] with their mean in (0, 1).	required
rate_b	float	Distinct probabilities on [0, 1] with their mean in (0, 1).	required
z_critical	float	Positive target magnitude of the normal statistic.	required

Returns

Name	Type	Description
	float	Unrounded count per group; callers round up when planning whole observations.

Notes

Equation: n = 2*p*(1-p)*(z_critical/(rate_a-rate_b))**2, p = (rate_a+rate_b)/2.

Assumptions: Equal-size independent Bernoulli groups and fixed hypothesized rates.

Excludes: Power beyond reaching a z threshold, rate-estimation uncertainty, and input validation.

References

Derivation: Rearrange z = gap/sqrt(2*p*(1-p)/n) for n.

Examples

>>> n = calc_equal_n_for_z(0.6, 0.4, 2.0)
>>> math.isclose(abs(calc_two_proportion_z(0.6, 0.4, n, n)), 2.0)
True

calc_expected_calibration_error

physics.statistics.calc_expected_calibration_error(bins)

Compute count-weighted absolute calibration error.

Expected calibration error of binned predictions: sum |B_i| / n * |acc_i - conf_i|.

bins is a sequence of (count, confidence, accuracy); the weight of each bin is its share of the predictions (Pakdaman Naeini et al. 2015; Guo et al. 2017, eq. 3).

Parameters

Name	Type	Description	Default
bins	sequence	Sequence of (count, confidence, accuracy) triples; rates in [0, 1].	required

Returns

Name	Type	Description
	float	Dimensionless weighted absolute probability error.

Notes

Equation: ECE = sum(count_i/total * abs(accuracy_i-confidence_i)).

Assumptions: Nonnegative bin counts with positive total; rates on the same probability scale.

Excludes: Bin selection, debiasing, and confidence intervals.

References

Derivation: Guo et al. (2017), On Calibration of Modern Neural Networks: count-weighted bin accuracy-confidence differences.

Examples

>>> calc_expected_calibration_error([(10, 0.5, 0.5)]) == 0
True

calc_exponential_decay_rate

physics.statistics.calc_exponential_decay_rate(initial, later, elapsed)

Fit an exponential rate through two readings.

Rate gamma of A(t) = A0 exp(-gamma t) through two readings: gamma = ln(A0 / A1) / t. elapsed may be a Quantity (the rate is then per its unit) or a plain number.

Parameters

Name	Type	Description	Default
initial	float	Positive initial reading, in the same units as later.	required
later	float	Positive later reading, in the same units as initial.	required
elapsed	float	Nonzero elapsed time, as a float or Pint quantity.	required

Returns

Name	Type	Description
	float or Quantity	Inverse-time rate, retaining the inverse units of elapsed when supplied.

Notes

Equation: gamma = log(initial/later)/elapsed.

Assumptions: Positive readings in matching units and nonzero elapsed time.

Excludes: Measurement noise and multiple-timescale decay.

References

Derivation: Take logarithms of later = initial*exp(-gamma*elapsed).

Examples

>>> calc_exponential_decay_rate(2, 2, 1) == 0
True

calc_fisher_exact_two_sided

physics.statistics.calc_fisher_exact_two_sided(a, b, c, d)

Compute a two-sided Fisher exact contingency-table p-value.

Parameters

Name	Type	Description	Default
a	int	Nonnegative cell counts of [[a, b], [c, d]], with nonzero total; callers must supply valid integer counts.	required
b	int	Nonnegative cell counts of [[a, b], [c, d]], with nonzero total; callers must supply valid integer counts.	required
c	int	Nonnegative cell counts of [[a, b], [c, d]], with nonzero total; callers must supply valid integer counts.	required
d	int	Nonnegative cell counts of [[a, b], [c, d]], with nonzero total; callers must supply valid integer counts.	required

Returns

Name	Type	Description
	float	Dimensionless probability of fixed-margin tables no more probable than the observed table.

Notes

Equation: sum_x P(X=x) for P(X=x) <= P(X=a)*(1+1e-9), with X hypergeometric.

Assumptions: Fixed row and column margins and independent observations under the null odds ratio of one.

Excludes: Alternative two-sided definitions, confidence intervals, and automatic input-domain checks.

References

Fisher, Statistical Methods for Research Workers (1934 edition), section 21.02. SciPy documents the same probability-ordering convention: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.fisher_exact.html

Examples

>>> math.isclose(calc_fisher_exact_two_sided(1, 1, 1, 1), 1.0, rel_tol=0.0, abs_tol=1e-12)
True

calc_gaussian_kl

physics.statistics.calc_gaussian_kl(mean_p, sd_p, mean_q, sd_q)

Compute the KL divergence between scalar normal distributions.

Parameters

Name	Type	Description	Default
mean_p	float	Mean of P, in the shared scalar unit.	required
sd_p	float	Positive standard deviation of P, in the shared scalar unit.	required
mean_q	float	Mean of Q, in the same unit as mean_p.	required
sd_q	float	Positive standard deviation of Q, in the shared scalar unit.	required

Returns

Name	Type	Description
	float	KL(P || Q), dimensionless and expressed in nats.

Notes

Equation: ln(sd_q/sd_p) + (sd_p**2 + (mean_p-mean_q)**2)/(2*sd_q**2) - 1/2.

Assumptions: The scalar inputs describe normal distributions, and logarithms are natural.

Excludes: Multivariate covariance and uncertainty in estimated distribution parameters.

References

Derivation: Substitute the two normal density functions into E_P[log(p/q)] and use E_P[(X-mean_q)**2] = sd_p**2 + (mean_p-mean_q)**2.

Examples

>>> calc_gaussian_kl(0, 1, 0, 1) == 0
True

calc_gaussian_kl_terms

physics.statistics.calc_gaussian_kl_terms(mean_p, sd_p, mean_q, sd_q)

Decompose the KL divergence between scalar normal distributions.

Parameters

Name	Type	Description	Default
mean_p	float	Mean of P, in the shared scalar unit.	required
sd_p	float	Positive standard deviation of P, in the shared scalar unit.	required
mean_q	float	Mean of Q, in the same unit as mean_p.	required
sd_q	float	Positive standard deviation of Q, in the shared scalar unit.	required

Returns

Name	Type	Description
	tuple of float	Log ratio (dimensionless), ratio numerator (squared input unit), and ratio denominator (squared input unit); the final -1/2 term is separate.

Notes

Equation: ln(sd_q/sd_p) + (sd_p**2 + (mean_p-mean_q)**2)/(2*sd_q**2) - 1/2.

Assumptions: P and Q are univariate normal distributions with positive standard deviations.

Excludes: Covariance matrices, estimation uncertainty, and nonnormal distribution shape.

References

Derivation: Substitute the two normal density functions into E_P[log(p/q)] and use E_P[(X-mean_q)**2] = sd_p**2 + (mean_p-mean_q)**2.

Examples

>>> log_term, numerator, denominator = calc_gaussian_kl_terms(0, 1, 0, 1)
>>> log_term + numerator / denominator - 0.5 == 0
True

calc_importance_sampling_exponential

physics.statistics.calc_importance_sampling_exponential(
    threshold,
    operational_mean,
    proposal_mean,
)

Compute exponential-tail importance-sampling moments.

Parameters

Name	Type	Description	Default
threshold	float or Quantity	Positive failure threshold gamma; same dimension as both means.	required
operational_mean	float or Quantity	Positive exponential mean m0, in threshold units.	required
proposal_mean	float or Quantity	Proposal mean mq > m0/2; plain floats must use a common unit.	required

Returns

Name	Type	Description
	dict of str to float	Dimensionless failure_probability, second_moment, relative_variance, ess_per_draw, and proposal_failure_share.

Notes

Equation: For lambda=1/m0 and eta=1/mq, P=exp(-lambda*gamma) and M2=lambda^2*exp(-(2*lambda-eta)*gamma)/(eta*(2*lambda-eta)); relative_variance=M2/P^2-1.

Assumptions: Exponential operational and proposal densities; likelihood-ratio weighting and finite second moment.

Excludes: Random draws, parameter fitting, proposal adaptation, and non-exponential tail behavior.

Operational model p = Exp(mean m0), failure when D > gamma, proposal q = Exp(mean m_q) with likelihood ratio w = p/q. Writing lam = 1/m0 and eta = 1/m_q (eta < 2 lam so the second moment is finite):

P_F                 = exp(-lam gamma)
E_q[w^2 1{D>gamma}] = lam^2 / (eta (2 lam - eta)) * exp(-(2 lam - eta) gamma)
relative variance   = E_q[w^2 1{D>gamma}] / P_F^2 - 1      (per draw)
ESS per draw        = E_q[w]^2 / E_q[w^2] = eta (2 lam - eta) / lam^2
P_q(D > gamma)      = exp(-eta gamma)

Draws for a relative standard error r are relative_variance / r**2. threshold, operational_mean, and proposal_mean are Pint quantities of one dimension (or plain floats in one unit).

Source: likelihood-ratio estimator, its second moment, and the effective sample size (sum w)^2 / sum w^2 from Owen 2013, ch. 9; the exponential closed forms follow by integrating those expressions.

References

Derivation: integrate (p(x)/q(x))^2*q(x) from gamma to infinity. Integrating from zero gives the weight second moment used for ess_per_draw.

Examples

>>> result = calc_importance_sampling_exponential(1.0, 1.0, 1.0)
>>> result["ess_per_draw"]
1.0

calc_jensen_shannon_divergence

physics.statistics.calc_jensen_shannon_divergence(p, q)

Measure divergence between aligned probability distributions.

Jensen-Shannon divergence between two aligned distributions, in nats.

JS = KL(p || m) / 2 + KL(q || m) / 2 with m = (p + q) / 2. It is symmetric and bounded by ln 2, which two distributions with disjoint support reach. Source: Lin, “Divergence Measures Based on the Shannon Entropy,” IEEE Trans. Inf. Theory 37(1), 1991.

Parameters

Name	Type	Description	Default
p	sequence	Aligned normalized nonnegative probabilities, dimensionless.	required
q	sequence	Aligned normalized nonnegative probabilities, dimensionless.	required

Returns

Name	Type	Description
	float	Divergence in nats.

Notes

Equation: JS = (KL(p||m)+KL(q||m))/2, m=(p+q)/2.

Assumptions: Aligned, nonnegative, normalized probabilities; natural logarithms.

Excludes: Automatic normalization, smoothing, and sampling uncertainty.

References

Derivation: Lin (1991), Divergence Measures Based on the Shannon Entropy; mixture-relative entropy definition.

Examples

>>> calc_jensen_shannon_divergence([0.5, 0.5], [0.5, 0.5]) == 0
True

calc_ks_coefficient

physics.statistics.calc_ks_coefficient(significance)

Compute an asymptotic two-sample KS critical coefficient.

Parameters

Name	Type	Description	Default
significance	float	Two-sided significance probability alpha in (0, 1).	required

Returns

Name	Type	Description
	float	Dimensionless large-sample coefficient for calc_ks_critical_value.

Notes

Equation: c(alpha) = sqrt(-log(alpha/2)/2).

Assumptions: The leading asymptotic KS tail approximation is adequate for the selected alpha.

Excludes: Exact finite-sample probabilities, ties, and parameter-estimation corrections.

About 1.36 at 0.05 (Ops.Monitoring.KsTestCoefficient) and 1.63 at 0.01; calc_ks_critical_value scales it by the two sample sizes.

References

Derivation: invert the leading tail approximation alpha=2*exp(-2*c**2).

Examples

>>> calc_ks_coefficient(0.01) > calc_ks_coefficient(0.05)
True

calc_ks_critical_value

physics.statistics.calc_ks_critical_value(coefficient, n_reference, n_current)

Scale a two-sample Kolmogorov-Smirnov critical coefficient.

Parameters

Name	Type	Description	Default
coefficient	float	Positive dimensionless asymptotic coefficient for the significance level.	required
n_reference	int	Positive independent reference-sample count.	required
n_current	int	Positive independent current-sample count.	required

Returns

Name	Type	Description
	float	Dimensionless critical maximum empirical-CDF difference.

Notes

Equation: D_crit = c*sqrt((n_reference+n_current)/(n_reference*n_current)).

Assumptions: Large independent samples from continuous distributions under the null.

Excludes: Exact finite-sample calibration, tied observations, and fitted-null corrections.

D_crit = c * sqrt((n_0 + n_t) / (n_0 * n_t)), where c is the coefficient for the chosen significance level (1.36 at 0.05, Ops.Monitoring.KsTestCoefficient).

References

Derivation: the asymptotic two-sample KS statistic multiplies D by sqrt(n*m/(n+m)); solve for D at the supplied critical coefficient.

Examples

>>> a = calc_ks_critical_value(1.0, 10, 10)
>>> b = calc_ks_critical_value(1.0, 40, 40)
>>> a == 2 * b
True

calc_ks_normal_location_shift

physics.statistics.calc_ks_normal_location_shift(shift_sd)

Compute the CDF distance for an equal-variance normal shift.

Parameters

Name	Type	Description	Default
shift_sd	float	Mean separation measured in shared standard-deviation units; either sign is allowed.	required

Returns

Name	Type	Description
	float	Dimensionless population KS distance in [0, 1].

Notes

Equation: D = erf(abs(shift_sd)/(2*sqrt(2))).

Assumptions: Two Gaussian populations share a positive variance.

Excludes: Empirical sampling noise, variance changes, and significance testing.

The CDF gap is largest midway between the means: 2 * Phi(shift_sd / 2) - 1, which equals erf(shift_sd / (2 * sqrt(2))).

References

Derivation: equal-density crossings occur at the midpoint of the means; evaluate the two normal CDFs there and use 2*Phi(z)-1=erf(z/sqrt(2)).

Examples

>>> calc_ks_normal_location_shift(0.0)
0.0

calc_linf_corner_l2

physics.statistics.calc_linf_corner_l2(linf_budget, coordinates)

Convert an infinity-norm corner radius to its Euclidean norm.

Parameters

Name	Type	Description	Default
linf_budget	float or Quantity	Nonnegative coordinate-wise perturbation bound, in the coordinate unit.	required
coordinates	int	Nonnegative number of equal-scale coordinates.	required

Returns

Name	Type	Description
	float or Quantity	Euclidean norm in the same unit as linf_budget.

Notes

Equation: l2 = linf_budget * sqrt(coordinates).

Assumptions: Each coordinate has absolute value equal to the bound; coordinate units agree.

Excludes: Typical perturbation size, weighted norms, and validation of the supplied domain.

References

Derivation: The sum of squared coordinates is coordinates * linf_budget**2; take its square root.

Examples

>>> calc_linf_corner_l2(0.0, 10) == 0
True

calc_mcnemar_exact_two_sided

physics.statistics.calc_mcnemar_exact_two_sided(first_only, second_only)

Compute the exact conditional McNemar p-value.

Exact two-sided McNemar p from the discordant counts (McNemar 1947).

Conditional on d = first_only + second_only discordant pairs, wins are Binomial(d, 1/2) under equal marginals: p = min(1, 2 sum_{j <= min(b, c)} C(d, j) 2^-d); 1 when d = 0.

Parameters

Name	Type	Description	Default
first_only	float	Nonnegative integer count accepted only by the first configuration.	required
second_only	float	Nonnegative integer count accepted only by the second configuration.	required

Returns

Name	Type	Description
	float	Dimensionless two-sided p-value.

Notes

Equation: p = min(1, 2*P[Binomial(b+c, 0.5) <= min(b,c)]).

Assumptions: Independent matched pairs and exchangeable discordant labels under the null.

Excludes: Effect-size intervals and dependence between pairs.

References

Derivation: Condition on discordance; equal marginals give a fair binomial count.

Examples

>>> calc_mcnemar_exact_two_sided(0, 0) == 1
True

calc_mcnemar_max_discordant

physics.statistics.calc_mcnemar_max_discordant(net_gap, critical)

Bound discordance for an uncorrected McNemar threshold.

Largest disagreement count b + c at which McNemar’s statistic still exceeds critical.

The statistic (no continuity correction) is (b - c)^2 / (b + c); with the net gap b - c fixed, it exceeds the 1-dof chi-square critical value only when b + c <= floor(gap^2 / critical) (McNemar 1947).

Parameters

Name	Type	Description	Default
net_gap	float	Signed integer difference between discordant cell counts.	required
critical	float	Positive dimensionless chi-square statistic threshold.	required

Returns

Name	Type	Description
	int	Maximum algebraic discordant count, before contingency-table feasibility.

Notes

Equation: d_max = floor(net_gap^2/critical).

Assumptions: Uncorrected chi-square approximation and a fixed net gap.

Excludes: Continuity correction, exact p-values, and parity feasibility of counts.

References

Derivation: Rearrange gap^2/d >= critical; equality is included by the returned floor.

Examples

>>> calc_mcnemar_max_discordant(4, 4) == 4
True

calc_mean_half_width

physics.statistics.calc_mean_half_width(sd, n_samples, confidence)

Calculate a Student t interval half-width for a mean.

Half-width of a confidence interval on a mean of n_samples values: t_{n-1} sd / sqrt(n).

Parameters

Name	Type	Description	Default
sd	float	Nonnegative standard deviation in the measurement unit.	required
n_samples	float	Observation count, at least two.	required
confidence	float	Central confidence probability in (0, 1).	required

Returns

Name	Type	Description
	float	Half-width in the same units as sd.

Notes

Equation: h = t_(n-1)*sd/sqrt(n_samples).

Assumptions: Independent observations; t coverage is exact for a normal population.

Excludes: Serial dependence and systematic measurement bias.

References

Derivation: Multiply the standard error by the two-sided Student t quantile.

Examples

>>> calc_mean_half_width(0, 10, 0.95) == 0
True

calc_monte_carlo_sample_size

physics.statistics.calc_monte_carlo_sample_size(
    failure_probability,
    relative_error,
)

Estimate independent draws for a relative standard-error target.

Parameters

Name	Type	Description	Default
failure_probability	float	Bernoulli failure probability p in (0, 1), dimensionless.	required
relative_error	float	Positive target relative standard error r, dimensionless.	required

Returns

Name	Type	Description
	float	Unrounded independent draw count, not a confidence-interval coverage guarantee.

Notes

Equation: N = (1-p)/(p*r^2).

Assumptions: Independent Bernoulli trials with known planning probability p.

Excludes: Confidence quantiles, correlated draws, and model misspecification.

Equation: N = (1 - p) / (p * r^2)

The binomial estimator p_hat = k/N has relative standard error sqrt((1 - p) / (N p)); solving for N at a target r gives the count.

Source: Owen, Monte Carlo Theory, Methods and Examples (2013): the variance of a sample mean of Bernoulli draws, as used for rare events in ch. 9.

References

Derivation: Var(p_hat) = p*(1-p)/N; divide its square root by p and solve for N. Owen (2013), Monte Carlo Theory, Methods and Examples, gives the sample-mean framework.

Examples

>>> calc_monte_carlo_sample_size(0.5, 0.5)
4.0

calc_multilegged_confidence

physics.statistics.calc_multilegged_confidence(leg_confidences)

Combine independent supporting legs by disjunction.

Parameters

Name	Type	Description	Default
leg_confidences	sequence of float	Nonempty probabilities p_i in [0, 1].	required

Returns

Name	Type	Description
	float	Dimensionless probability that at least one sufficient leg holds.

Notes

Equation: P(any leg) = 1-product(1-p_i).

Assumptions: Each leg suffices alone and their failures are independent.

Excludes: Common-cause dependence and evidential calibration of the input probabilities.

Equation: P(claim) = 1 - prod_i (1 - p_i)

Holds only if the legs fail independently; common-cause dependence requires a joint probability model instead of this product.

Source: Bloomfield and Littlewood 2003 (DSN), multi-legged arguments under the independence assumption.

References

Derivation: all legs fail with the product of their failure probabilities; take its complement.

Examples

>>> calc_multilegged_confidence([0.5, 0.5])
0.75

calc_newcombe_paired_difference

physics.statistics.calc_newcombe_paired_difference(
    both,
    first_only,
    second_only,
    neither,
    z_critical,
)

Estimate a paired rate difference and its interval.

Paired difference of two acceptance rates with Newcombe’s method-10 interval (Newcombe 1998).

The four cells count the tasks both configurations accept, only the first, only the second, and neither. Returns (delta, low, high) for the first configuration’s rate minus the second’s. Each marginal Wilson interval is combined through the continuity-corrected phi of the 2x2 table, so positive within-pair dependence narrows the interval.

Parameters

Name	Type	Description	Default
both	float	Nonnegative integer count of jointly accepted pairs.	required
first_only	float	Nonnegative integer count accepted only by the first configuration.	required
second_only	float	Nonnegative integer count accepted only by the second configuration.	required
neither	float	Nonnegative integer count accepted by neither configuration.	required
z_critical	float	Positive dimensionless normal quantile.	required

Returns

Name	Type	Description
	tuple of float	Difference, lower limit, upper limit on the dimensionless rate scale.

Notes

Equation: delta = p1-p2; marginal Wilson widths combine with corrected phi covariance.

Assumptions: Independent pairs and nonnegative integer contingency counts.

Excludes: Unpaired sampling and repeated-measures cluster corrections.

References

Derivation: Newcombe (1998), method 10: combine marginal Wilson score limits using the paired-table correlation.

Examples

>>> calc_newcombe_paired_difference(5, 2, 2, 5, 1.96)[0] == 0
True

calc_noncentral_gate_detection

physics.statistics.calc_noncentral_gate_detection(gate, dof, noncentrality)

Compute the probability that a biased statistic crosses a gate.

Parameters

Name	Type	Description	Default
gate	float	Positive dimensionless squared-statistic threshold.	required
dof	int	Positive number of degrees of freedom.	required
noncentrality	float	Nonnegative dimensionless squared normalized bias; b**2/S in the scalar case.	required

Returns

Name	Type	Description
	float	Probability in [0, 1] of exceeding the gate.

Notes

Equation: detection = 1 - F_ncx2(gate; dof, noncentrality), evaluated by the survival function.

Assumptions: Gaussian residuals with fixed covariance and fixed normalized bias.

Excludes: Time-varying bias, multiple-update detection time, and uncertainty in covariance or bias.

References

Bar-Shalom, Li, and Kirubarajan (2001), sections 1.4.17 and 5.4. https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.ncx2.html

Examples

>>> calc_noncentral_gate_detection(5.0, 2, 4.0) > calc_noncentral_gate_detection(5.0, 2, 0.0)
True

calc_one_proportion_sample_size

physics.statistics.calc_one_proportion_sample_size(
    null_rate,
    alternative_rate,
    z_alpha,
    z_beta,
)

Plan a one-sample binomial comparison of two rates.

Outcomes a one-sample binomial test needs to tell alternative_rate from null_rate.

Normal approximation with the variance of each hypothesis on its own side: n = ((z_alpha sqrt(p0 (1 - p0)) + z_beta sqrt(p1 (1 - p1))) / (p0 - p1))^2. The count grows as one over the squared gap, so halving the gap needs about four times the outcomes.

Parameters

Name	Type	Description	Default
null_rate	float	Null Bernoulli probability in (0, 1).	required
alternative_rate	float	Distinct alternative Bernoulli probability in (0, 1).	required
z_alpha	float	Dimensionless normal significance quantile; default 1.96 where specified.	required
z_beta	float	Dimensionless normal power quantile; default 0.84 where specified.	required

Returns

Name	Type	Description
	float	Unrounded count; round upward for planning.

Notes

Equation: n = ((z_alpha*sqrt(p0*(1-p0))+z_beta*sqrt(p1*(1-p1)))/(p0-p1))^2.

Assumptions: Independent Bernoulli observations and distinct interior hypothesis rates.

Excludes: Exact binomial power and sequential sampling.

References

Derivation: Balance normal null and alternative standard errors against the rate gap.

Examples

>>> calc_one_proportion_sample_size(0.5, 0.6, 2, 1) > 0
True

calc_optimal_exponential_proposal_mean

physics.statistics.calc_optimal_exponential_proposal_mean(
    threshold,
    operational_mean,
)

Minimize the exponential importance-sampling second moment.

Parameters

Name	Type	Description	Default
threshold	float or Quantity	Positive tail threshold gamma.	required
operational_mean	float or Quantity	Positive operational exponential mean, in threshold units.	required

Returns

Name	Type	Description
	float or Quantity	Optimal exponential proposal mean in threshold units.

Notes

Equation: With lambda=1/m0, solve gamma*eta^2-(2+2*lambda*gamma)*eta+2*lambda=0 for the smaller root, then return 1/eta.

Assumptions: The proposal family is exponential and minimizes the same tail-estimator second moment as calc_importance_sampling_exponential.

Excludes: Optimization over other proposal families and finite-sample adaptation costs.

Minimizing the second moment of calc_importance_sampling_exponential over eta gives the root below 2 lam of

gamma eta^2 - (2 + 2 lam gamma) eta + 2 lam = 0,

and the proposal mean is 1/eta. Returned in the units of threshold.

Source: first-order condition on the exponential second moment (Owen 2013, ch. 9, exponential tilting of a light-tailed tail).

References

Derivation: differentiate the logarithm of the exponential second moment with respect to proposal rate eta and solve the resulting quadratic.

Examples

>>> proposal = calc_optimal_exponential_proposal_mean(1.0, 1.0)
>>> proposal > 1.0
True

calc_paired_difference_variance

physics.statistics.calc_paired_difference_variance(
    discordant_fraction,
    difference,
)

Calculate variance of a signed paired outcome.

Per-pair variance of a paired candidate-minus-baseline outcome in {-1, 0, +1}: d - delta^2.

Parameters

Name	Type	Description	Default
discordant_fraction	float	Fraction of pairs with different outcomes, in [0, 1].	required
difference	float	Mean difference on the outcome scale; nonzero for sample planning.	required

Returns

Name	Type	Description
	float	Dimensionless per-pair variance.

Notes

Equation: V = discordant_fraction - difference^2.

Assumptions: Outcomes are -1, 0, or +1 with compatible supplied moments.

Excludes: Variance estimation uncertainty and clustering.

References

Derivation: For signed discordance X, E[X^2] is the discordant fraction.

Examples

>>> calc_paired_difference_variance(0.5, 0) == 0.5
True

calc_poisson_exact_interval

physics.statistics.calc_poisson_exact_interval(events, alpha=0.05)

Invert Poisson tails to bound an event-count mean.

Parameters

Name	Type	Description	Default
events	int	Nonnegative observed event count; implementation coerces with int, so callers should supply integers.	required
alpha	float	Two-sided tail error probability in (0, 1); default 0.05.	0.05

Returns

Name	Type	Description
	tuple of float	Lower and upper Poisson mean counts; divide by exposure to obtain rate units.

Notes

Equation: P_mu(N >= k) = alpha/2 at lower; P_mu(N <= k) = alpha/2 at upper.

Assumptions: A homogeneous Poisson count and sufficiently moderate counts for direct probability recurrence.

Excludes: Overdispersion, exposure uncertainty, extreme-count numerical stabilization, and validation of count or alpha.

References

Garwood (1936), Fiducial Limits for the Poisson Distribution, Biometrika 28(3/4). Derivation: Invert the equal-tail tests of the Poisson mean; the zero-count lower endpoint is zero.

Examples

>>> lower, upper = calc_poisson_exact_interval(0)
>>> lower == 0 and upper > 0
True

calc_poisson_mean_for_reserve

physics.statistics.calc_poisson_mean_for_reserve(reserve, target, iterations=60)

Find a mean-count boundary for a fixed Poisson reserve.

Parameters

Name	Type	Description	Default
reserve	int	Positive integer reserve count.	required
target	float	Upper-tail probability target in (0, 1).	required
iterations	int	Positive bisection iteration count; default 60.	60

Returns

Name	Type	Description
	float	Lower-side approximation to the mean with tail probability equal to target.

Notes

Equation: Solve P(Poisson(mean)>=reserve)=target by expanding the upper bracket then bisecting.

Assumptions: The Poisson upper tail increases monotonically with its mean; valid reserve and target guarantee a crossing.

Excludes: Exact root arithmetic, extreme-tail numerical stabilization, and non-Poisson demand.

Bisection on the tail, which rises monotonically with the mean.

References

Derivation: monotonic inversion of the Poisson survival probability with bisection.

Examples

>>> mean = calc_poisson_mean_for_reserve(1, 0.5)
>>> abs(calc_poisson_tail(mean, 1) - 0.5) < 1e-12
True

calc_poisson_reserve

physics.statistics.calc_poisson_reserve(mean, target)

Find the smallest reserve satisfying a Poisson tail target.

Parameters

Name	Type	Description	Default
mean	float	Nonnegative expected count during replenishment lead time.	required
target	float	Required upper-tail probability ceiling in (0, 1).	required

Returns

Name	Type	Description
	int	Smallest reserve s with P(N >= s) < target; equality does not pass.

Notes

Equation: s = min{integer s>=0 : P(Poisson(mean)>=s) < target}.

Assumptions: Poisson arrivals and fixed replenishment lead time, with a valid positive target.

Excludes: Variable lead times, correlated failures, and a stable extreme-tail evaluation.

References

Derivation: enumerate integer reserves in increasing order and stop at the first upper tail below the target.

Examples

>>> calc_poisson_reserve(0.0, 0.1)
1

calc_poisson_tail

physics.statistics.calc_poisson_tail(mean, at_least)

Compute the probability of at least a specified event count.

Parameters

Name	Type	Description	Default
mean	float	Nonnegative expected Poisson count in one observation interval.	required
at_least	int	Integer event-count threshold; nonpositive thresholds return one.	required

Returns

Name	Type	Description
	float	Dimensionless upper-tail probability, evaluated by direct CDF summation.

Notes

Equation: P(N>=s) = 1-sum(exp(-mean)*mean**k/k!, k=0,…,s-1).

Assumptions: Poisson counts in a fixed exposure interval.

Excludes: Overdispersion, nonstationary intensity, and cancellation-resistant extreme-tail numerics.

The chance that at least s events fall in an interval whose expected count is mean; 1 for s <= 0. A pool with s standby instances and a constant replenishment lead time misses a request exactly when at least s earlier replacements are still outstanding.

References

Derivation: complement the sum of the Poisson mass function below s.

Examples

>>> calc_poisson_tail(0.0, 1)
0.0

calc_population_stability_index

physics.statistics.calc_population_stability_index(
    expected,
    actual,
    epsilon=1e-12,
)

Measure drift between two aligned binned distributions.

Parameters

Name	Type	Description	Default
expected	sequence of float	Aligned, equal-length nonnegative bin probabilities, each summing to one.	required
actual	sequence of float	Aligned, equal-length nonnegative bin probabilities, each summing to one.	required
epsilon	float	Positive dimensionless probability floor; default 1e-12 avoids log of zero.	1e-12

Returns

Name	Type	Description
	float	Dimensionless symmetrized divergence after flooring; zero for identical bins.

Notes

Equation: PSI = sum((a_i - e_i) * log(a_i / e_i)), after flooring each probability.

Assumptions: Both distributions use identical bins; flooring does not renormalize them.

Excludes: Bin selection, uncertainty intervals, and automatic drift-threshold decisions.

Measures distribution drift between a reference (“expected”) and an observed (“actual”) binned distribution: PSI = sum_i (actual_i - expected_i) * ln(actual_i / expected_i) (the symmetrized KL divergence). Interpretation thresholds live in Ops.Monitoring (the single source of truth the book follows: 0.10 warn / 0.20 review / 0.25 critical) — do not restate bands here.

References

Derivation: adding KL(actual || expected) and KL(expected || actual) yields the stated sum; flooring is a numerical convention.

Examples

>>> calc_population_stability_index([0.5, 0.5], [0.5, 0.5])
0.0

calc_proportion_gap_interval

physics.statistics.calc_proportion_gap_interval(
    rate_a,
    rate_b,
    n_a,
    n_b,
    z_critical,
)

Bound a proportion difference using its unpooled standard error.

Parameters

Name	Type	Description	Default
rate_a	float	Observed independent-group proportions on [0, 1].	required
rate_b	float	Observed independent-group proportions on [0, 1].	required
n_a	int	Positive observation counts for the two groups.	required
n_b	int	Positive observation counts for the two groups.	required
z_critical	float	Positive normal critical value for the chosen coverage.	required

Returns

Name	Type	Description
	tuple of float	Lower and upper endpoints in probability-difference units, without clipping.

Notes

Equation: interval = (rate_a-rate_b) +/- z_critical*sqrt(rate_a*(1-rate_a)/n_a + rate_b*(1-rate_b)/n_b).

Assumptions: Independent groups and adequate sample counts for a normal approximation.

Excludes: Paired observations, exact small-count coverage, and restricting endpoints to the feasible interval.

References

Derivation: Add independent proportion variances, take the square root, and form a symmetric normal interval around the observed gap.

Examples

>>> lo, hi = calc_proportion_gap_interval(0.5, 0.5, 100, 100, 1.96)
>>> lo < 0 < hi
True

calc_proportion_gap_se

physics.statistics.calc_proportion_gap_se(rate_a, rate_b, n_a, n_b)

Estimate the standard error of two independent proportions.

Parameters

Name	Type	Description	Default
rate_a	float	Observed independent-group proportions on [0, 1].	required
rate_b	float	Observed independent-group proportions on [0, 1].	required
n_a	int	Positive observation counts for the two groups.	required
n_b	int	Positive observation counts for the two groups.	required

Returns

Name	Type	Description
	float	Nonnegative standard error in probability-difference units.

Notes

Equation: se = sqrt(rate_a*(1-rate_a)/n_a + rate_b*(1-rate_b)/n_b).

Assumptions: Each group is an independent Bernoulli sample; plug-in observed rates estimate its variance.

Excludes: Paired covariance, clustering, finite-population corrections, and input validation.

References

Derivation: Variances of independent sample means add when subtracting the means.

Examples

>>> calc_proportion_gap_se(0.5, 0.5, 400, 400) < calc_proportion_gap_se(0.5, 0.5, 100, 100)
True

calc_proportion_half_width

physics.statistics.calc_proportion_half_width(rate, n_samples, z_critical)

Estimate the normal interval half-width for a proportion.

Parameters

Name	Type	Description	Default
rate	float	Bernoulli probability estimate p on [0, 1].	required
n_samples	int	Positive number of independent samples.	required
z_critical	float	Nonnegative normal critical value for the desired coverage.	required

Returns

Name	Type	Description
	float	Nonnegative half-width in probability units.

Notes

Equation: half_width = z_critical*sqrt(rate*(1-rate)/n_samples).

Assumptions: Independent Bernoulli observations and a normal plug-in approximation.

Excludes: Boundary corrections, exact coverage, dependence, and input validation.

References

Derivation: Multiply the standard deviation of a Bernoulli sample mean by the chosen normal critical value.

Examples

>>> calc_proportion_half_width(0.5, 400, 1.96) < calc_proportion_half_width(0.5, 100, 1.96)
True

calc_sample_collection_time

physics.statistics.calc_sample_collection_time(
    samples,
    request_rate,
    traffic_share,
)

Estimate collection time for a fixed traffic allocation.

Parameters

Name	Type	Description	Default
samples	int or Quantity	Nonnegative target observation count.	required
request_rate	float or Quantity	Positive requests per time unit; plain floats set the returned time unit.	required
traffic_share	float	Allocated fraction in (0, 1].	required

Returns

Name	Type	Description
	float or Quantity	Time to collect samples at a constant effective observation rate.

Notes

Equation: t = samples/(request_rate*traffic_share).

Assumptions: Every allocated request supplies one observation and the traffic rate is constant.

Excludes: Burstiness, dropped observations, dependence, and ramp-up scheduling.

A canary stage that receives traffic_share of request_rate gathers its planned guardrail sample in t = n / (r p).

References

Derivation: observation count equals effective rate times elapsed time; solve for elapsed time.

Examples

>>> calc_sample_collection_time(100, 20.0, 0.5)
10.0

calc_samples_for_half_width

physics.statistics.calc_samples_for_half_width(rate, half_width, z_critical)

Plan independent Bernoulli draws for a target half-width.

Independent draws for a normal-approximation half-width at rate: ceil(z^2 p (1 - p) / h^2).

Parameters

Name	Type	Description	Default
rate	float	Bernoulli probability in [0, 1], dimensionless.	required
half_width	float	Positive absolute rate half-width, dimensionless.	required
z_critical	float	Positive dimensionless normal quantile.	required

Returns

Name	Type	Description
	int	Draw count, rounded upward.

Notes

Equation: n = ceil(z^2 p(1-p)/h^2).

Assumptions: Independent Bernoulli observations; normal approximation at the supplied rate.

Excludes: Exact binomial coverage and clustered observations.

References

Derivation: Invert the normal half-width z sqrt(p(1-p)/n).

Examples

>>> calc_samples_for_half_width(0.5, 0.1, 2) > 0
True

calc_samples_per_arm_for_difference

physics.statistics.calc_samples_per_arm_for_difference(
    sd,
    difference,
    z_alpha,
    z_beta,
)

Plan equal-arm observations for a mean difference.

Samples per arm for a two-sided two-sample test of a mean difference: 2 (z_a + z_b)^2 sd^2 / d^2.

The normal-approximation power formula; z_alpha is the two-sided critical value (1.96 at 95 percent) and z_beta the power quantile (0.84 at 80 percent power). Not rounded up.

Parameters

Name	Type	Description	Default
sd	float	Nonnegative standard deviation in the measurement unit.	required
difference	float	Mean difference on the outcome scale; nonzero for sample planning.	required
z_alpha	float	Dimensionless normal significance quantile; default 1.96 where specified.	required
z_beta	float	Dimensionless normal power quantile; default 0.84 where specified.	required

Returns

Name	Type	Description
	float	Unrounded per-arm sample count; round upward for planning.

Notes

Equation: n = 2*(z_alpha+z_beta)^2*sd^2/difference^2.

Assumptions: Independent arms with equal variance and a normal power approximation.

Excludes: Small-sample t iteration and unequal allocation.

References

Derivation: Solve the normal difference-of-means power equation for equal arm size.

Examples

>>> calc_samples_per_arm_for_difference(1, 1, 2, 1) > 0
True

calc_smoothing_radius

physics.statistics.calc_smoothing_radius(noise_sd, top_class_probability)

Compute the binary Gaussian randomized-smoothing radius.

Parameters

Name	Type	Description	Default
noise_sd	float or Quantity	Nonnegative Gaussian noise standard deviation, in the input norm unit; the function does not validate its sign.	required
top_class_probability	float	Certified lower probability bound for the top class, between 0.5 and one.	required

Returns

Name	Type	Description
	float or Quantity	Certified l2 radius in the same unit and scalar representation as noise_sd.

Notes

Equation: R = noise_sd * Phi**(-1)(top_class_probability).

Assumptions: Isotropic Gaussian smoothing and a valid probability lower bound; the binary complement bounds the competing class.

Excludes: Estimation of the probability bound, nonspherical noise, and certification of an unsmoothed classifier.

References

Cohen, Rosenfeld, and Kolter (2019), Certified Adversarial Robustness via Randomized Smoothing. https://arxiv.org/abs/1902.02918

Examples

>>> calc_smoothing_radius(2.0, 0.75) == 2 * calc_smoothing_radius(1.0, 0.75)
True

calc_spearman_rank_correlation

physics.statistics.calc_spearman_rank_correlation(xs, ys)

Correlate paired observations through their average ranks.

Spearman’s rank correlation: Pearson correlation of the average ranks (ties share a rank).

Parameters

Name	Type	Description	Default
xs	sequence	First nonconstant sequence of numeric observations.	required
ys	sequence	Second nonconstant sequence paired with xs.	required

Returns

Name	Type	Description
	float	Dimensionless rank correlation between minus one and one.

Notes

Equation: rho = centered rank cross-product / product of centered rank norms.

Assumptions: Equal-length paired sequences with at least two values and nonconstant ranks.

Excludes: P-values, missing-value handling, and causal interpretation.

References

Derivation: Replace values by average ranks, then apply the Pearson correlation formula.

Examples

>>> calc_spearman_rank_correlation([1, 2, 3], [2, 4, 6]) == 1
True

calc_split_conformal_min_count

physics.statistics.calc_split_conformal_min_count(alpha)

Find the minimum calibration count for a finite conformal rank.

Parameters

Name	Type	Description	Default
alpha	float	Desired miscoverage probability in (0, 1).	required

Returns

Name	Type	Description
	int	Minimum calibration-score count for which the one-based rank does not exceed the count.

Notes

Equation: n_min = ceil((1-alpha)/alpha).

Assumptions: The split-conformal rank is ceil((n+1)*(1-alpha)); arithmetic is evaluated in floating point.

Excludes: Conditional coverage, score construction, and guarantees under nonexchangeable data.

References

Derivation: Solve (n+1)*(1-alpha) <= n for n, then round upward. Angelopoulos and Bates: https://arxiv.org/abs/2107.07511

Examples

>>> n = calc_split_conformal_min_count(0.1)
>>> calc_split_conformal_rank(n, 0.1) <= n
True

calc_split_conformal_rank

physics.statistics.calc_split_conformal_rank(calibration_count, alpha)

Find the finite-sample split-conformal calibration rank.

Parameters

Name	Type	Description	Default
calibration_count	int	Number n of calibration scores; at least one.	required
alpha	float	Desired marginal miscoverage probability, between zero and one.	required

Returns

Name	Type	Description
	int	One-based rank; a value greater than n means no finite calibration quantile is available.

Notes

Equation: rank = ceil((n + 1) * (1 - alpha)).

Assumptions: Calibration scores and the new score are exchangeable; the scoring rule is fixed independently of calibration.

Excludes: Conditional coverage, distribution shift, and constructing the prediction set itself.

References

Angelopoulos and Bates, A Gentle Introduction to Conformal Prediction and Distribution-Free Uncertainty Quantification. https://arxiv.org/abs/2107.07511

Examples

>>> calc_split_conformal_rank(100, 0.1) <= 100
True

calc_sprt_boundaries

physics.statistics.calc_sprt_boundaries(alpha, beta)

Compute Wald stopping boundaries for a log-likelihood ratio.

Parameters

Name	Type	Description	Default
alpha	float	Allowed false acceptance of H1 under H0, between zero and one.	required
beta	float	Allowed false rejection of H1 under H1, between zero and one.	required

Returns

Name	Type	Description
	tuple of float	Upper acceptance boundary for H1, then lower rejection boundary, in log-likelihood nats.

Notes

Equation: upper = ln((1-beta)/alpha); lower = ln(beta/(1-alpha)).

Assumptions: Wald simple-hypothesis sequential test; alpha + beta < 1 gives the usual separated boundaries.

Excludes: Discrete overshoot correction and validation that alpha + beta is below one.

References

Wald (1945), Sequential Tests of Statistical Hypotheses. Derivation: Take logarithms of the likelihood-ratio approximations (1-beta)/alpha and beta/(1-alpha).

Examples

>>> upper, lower = calc_sprt_boundaries(0.05, 0.05)
>>> upper > 0 > lower
True

calc_sprt_expected_samples

physics.statistics.calc_sprt_expected_samples(
    upper,
    lower,
    drift,
    upper_probability,
)

Approximate a likelihood walk stopping count.

Wald’s expected stopping sample count, neglecting overshoot.

The walk ends at upper with probability upper_probability and at lower otherwise; dividing the expected terminal log-likelihood by the expected step drift gives E[N] (Wald 1945, average sample number).

Parameters

Name	Type	Description	Default
upper	float	Upper terminal log-likelihood boundary, dimensionless.	required
lower	float	Lower terminal log-likelihood boundary, dimensionless.	required
drift	float	Nonzero mean log-likelihood increment per observation.	required
upper_probability	float	Probability of terminating at the upper boundary, in [0, 1].	required

Returns

Name	Type	Description
	float	Approximate sample count; inputs must describe a consistent walk.

Notes

Equation: E[N] = (upper_probability*upper+(1-upper_probability)*lower)/drift.

Assumptions: A stopped walk satisfying Wald expectation conditions with nonzero mean increment.

Excludes: Boundary overshoot and finite-horizon truncation.

References

Derivation: Wald expectation identity equates expected terminal sum to drift times expected count.

Examples

>>> calc_sprt_expected_samples(2, -2, 1, 1) == 2
True

calc_sprt_increments

physics.statistics.calc_sprt_increments(p_unacceptable, p_acceptable)

Compute Bernoulli SPRT log-likelihood increments.

Parameters

Name	Type	Description	Default
p_unacceptable	float	Failure probability under H0, greater than p_acceptable and below one.	required
p_acceptable	float	Failure probability under H1, greater than zero and below p_unacceptable.	required

Returns

Name	Type	Description
	tuple of float	Success increment then failure increment, in log-likelihood nats.

Notes

Equation: success = ln((1-p1)/(1-p0)); failure = ln(p1/p0).

Assumptions: Independent Bernoulli trials and fixed simple hypotheses p0 and p1.

Excludes: Adaptive hypotheses, dependent outcomes, and stopping-rule overshoot corrections.

References

Wald (1945), Sequential Tests of Statistical Hypotheses. Derivation: Divide each Bernoulli event probability under H1 by its probability under H0, then take the natural logarithm.

Examples

>>> success, failure = calc_sprt_increments(0.2, 0.1)
>>> success > 0 > failure
True

calc_standard_error

physics.statistics.calc_standard_error(sd, n_samples)

Compute the standard error of an independent sample mean.

Parameters

Name	Type	Description	Default
sd	float or Quantity	Nonnegative population or estimated standard deviation, in observation units.	required
n_samples	int	Positive independent observation count.	required

Returns

Name	Type	Description
	float or Quantity	Mean standard error in the same units as sd.

Notes

Equation: SE = sd/sqrt(n_samples).

Assumptions: Independent equal-variance observations; sd describes their spread.

Excludes: Autocorrelation, unequal weights, confidence quantiles, and uncertainty in sd.

References

Derivation: the variance of an average of n independent equal-variance observations is sd**2/n.

Examples

>>> calc_standard_error(2.0, 4)
1.0

calc_t_critical

physics.statistics.calc_t_critical(degrees_of_freedom, confidence)

Return a two-sided Student t critical value.

Two-sided Student t critical value: the (1 + confidence) / 2 quantile (2.78 at 4 dof, 0.95).

Non-integer degrees of freedom (Welch-Satterthwaite) are allowed.

Parameters

Name	Type	Description	Default
degrees_of_freedom	float	Positive degrees of freedom, possibly fractional.	required
confidence	float	Central confidence probability in (0, 1).	required

Returns

Name	Type	Description
	float	Dimensionless positive t quantile.

Notes

Equation: t = inverse_t_cdf((1+confidence)/2, degrees_of_freedom).

Assumptions: Positive degrees of freedom; confidence lies between zero and one.

Excludes: Automatic choice of degrees of freedom from data.

References

Derivation: Invert the symmetric Student t distribution at the requested central coverage.

Examples

>>> calc_t_critical(10, 0.95) > 0
True

calc_task_averaged_variance

physics.statistics.calc_task_averaged_variance(
    between_task_variance,
    within_task_variance,
    tasks,
    runs_per_task,
)

Calculate variance of an equally weighted task average.

Variance of a task-averaged acceptance estimate: sigma_task^2 / N + sigma_run^2 / (N r).

Tasks are sampled independently and runs are independent given the task; more runs shrink only the second term, more tasks shrink both.

Parameters

Name	Type	Description	Default
between_task_variance	float	Nonnegative variance of task probabilities, dimensionless.	required
within_task_variance	float	Nonnegative conditional outcome variance, dimensionless.	required
tasks	float	Positive number of independently sampled tasks.	required
runs_per_task	float	Positive equal run count per task.	required

Returns

Name	Type	Description
	float	Dimensionless variance of the estimated average rate.

Notes

Equation: V = between_task_variance/tasks + within_task_variance/(tasks*runs_per_task).

Assumptions: Independent tasks and conditionally independent runs with equal replication.

Excludes: Correlated runs and unequal task weights.

References

Derivation: Apply total variance to the mean of replicated task outcomes.

Examples

>>> calc_task_averaged_variance(1, 1, 10, 2) < calc_task_averaged_variance(1, 1, 10, 1)
True

calc_two_proportion_detectable_lift

physics.statistics.calc_two_proportion_detectable_lift(
    baseline_rate,
    n_per_variant,
    z_alpha=1.96,
    z_beta=0.84,
)

Invert the common-variance two-arm planning formula.

Smallest absolute lift a two-arm test with n_per_variant per arm can detect.

The common-variance sample-size formula of calc_two_proportion_sample_size solved for the effect: delta = sqrt(2 (z_alpha + z_beta)^2 p (1 - p) / n).

Parameters

Name	Type	Description	Default
baseline_rate	float	Baseline probability in (0, 1).	required
n_per_variant	float	Positive observation count per arm.	required
z_alpha	float	Dimensionless normal significance quantile; default 1.96 where specified.	1.96
z_beta	float	Dimensionless normal power quantile; default 0.84 where specified.	0.84

Returns

Name	Type	Description
	float	Absolute rate difference, dimensionless.

Notes

Equation: delta = sqrt(2*(z_alpha+z_beta)^2*p*(1-p)/n_per_variant).

Assumptions: Independent equally sized arms and baseline variance reused for both.

Excludes: Exact power, unequal allocation, and pooled alternative variance.

References

Derivation: Solve the normal sample-size equation for its absolute effect.

Examples

>>> calc_two_proportion_detectable_lift(0.5, 200) < calc_two_proportion_detectable_lift(0.5, 100)
True

calc_two_proportion_sample_size

physics.statistics.calc_two_proportion_sample_size(
    baseline_rate,
    detectable_lift,
    z_alpha=1.96,
    z_beta=0.84,
    pooled=False,
)

Estimate per-arm sample size for a two-proportion comparison.

Parameters

Name	Type	Description	Default
baseline_rate	float	Control probability p in (0, 1).	required
detectable_lift	float	Positive absolute probability difference delta; pooled mode also requires p + delta < 1.	required
z_alpha	float	Positive normal significance quantile; default 1.96 is the usual rounded two-sided 5-percent value.	1.96
z_beta	float	Normal power quantile; default 0.84 is the usual rounded 80-percent value.	0.84
pooled	bool	Use both arm variances when True; default False reuses baseline variance.	False

Returns

Name	Type	Description
	float	Unrounded sample count per arm; round upward before planning collection.

Notes

Equation: n = 2*(z_alpha+z_beta)^2*p*(1-p)/delta^2 by default. Pooled mode uses (z_alpha*sqrt(2*p_bar*(1-p_bar)) + z_beta*sqrt(p*(1-p)+(p+delta)*(1-p-delta)))^2/delta^2.

Assumptions: Independent Bernoulli outcomes, equal allocation, and a normal approximation.

Excludes: Sequential monitoring, multiple-testing adjustment, clustering, and exact small-sample power.

Standard normal-approximation power formula: n = 2 * (z_alpha + z_beta)^2 * p(1-p) / delta^2, using the baseline rate’s variance for both arms (a common simplification valid for small lifts). Defaults correspond to a two-sided 5% significance level (z = 1.96) and 80% power (z = 0.84).

References

Derivation: a normal critical separation plus the desired power separation is equated to the mean difference; solving the variance-of-two-means expression gives n.

Examples

>>> small = calc_two_proportion_sample_size(0.5, 0.02)
>>> large = calc_two_proportion_sample_size(0.5, 0.01)
>>> large == small * 4
True

calc_two_proportion_z

physics.statistics.calc_two_proportion_z(rate_a, rate_b, n_a, n_b)

Compute a pooled two-proportion z statistic.

Parameters

Name	Type	Description	Default
rate_a	float	Observed independent-group proportions on [0, 1].	required
rate_b	float	Observed independent-group proportions on [0, 1].	required
n_a	int	Positive observation counts for the two groups.	required
n_b	int	Positive observation counts for the two groups.	required

Returns

Name	Type	Description
	float	Signed dimensionless normal statistic for rate_a minus rate_b.

Notes

Equation: p = (n_a*rate_a+n_b*rate_b)/(n_a+n_b); z = (rate_a-rate_b)/sqrt(p*(1-p)*(1/n_a+1/n_b)).

Assumptions: Independent Bernoulli groups, a common rate under the null, and pooled p in (0, 1).

Excludes: Small-sample exact inference, continuity correction, dependence, and input validation.

References

Derivation: Estimate the common null probability by the pooled count and standardize the difference using the sum of independent Bernoulli-mean variances.

Examples

>>> calc_two_proportion_z(0.5, 0.5, 100, 100) == 0
True

calc_wald_interval

physics.statistics.calc_wald_interval(successes, trials, z_critical)

Compute an unclipped Wald interval for a binomial proportion.

Parameters

Name	Type	Description	Default
successes	int	Accepted outcomes x, from zero through trials.	required
trials	int	Positive number n of independent Bernoulli trials.	required
z_critical	float	Nonnegative normal critical value; callers choose it for the desired coverage.	required

Returns

Name	Type	Description
	tuple of float	Lower and upper probability endpoints; they may fall outside [0, 1].

Notes

Equation: p_hat +/- z_critical*sqrt(p_hat*(1-p_hat)/trials), p_hat = successes/trials.

Assumptions: Independent identically distributed Bernoulli trials and a normal plug-in approximation.

Excludes: Boundary corrections and reliable small-count coverage; zero or all successes produce zero width.

References

Brown, Cai, and DasGupta (2001), Interval Estimation for a Binomial Proportion. Derivation: Use the observed Bernoulli variance in the symmetric normal interval for its sample mean.

Examples

>>> calc_wald_interval(0, 10, 1.96) == (0.0, 0.0)
True

calc_welch_interval

physics.statistics.calc_welch_interval(
    mean_a,
    sd_a,
    n_a,
    mean_b,
    sd_b,
    n_b,
    confidence,
)

Construct a Welch interval for two independent means.

Welch’s interval on mean_a - mean_b for two samples with unequal spreads (Welch 1947).

se = sqrt(s_a^2 / n_a + s_b^2 / n_b); the Welch-Satterthwaite degrees of freedom are se^4 / ((s_a^2/n_a)^2 / (n_a - 1) + (s_b^2/n_b)^2 / (n_b - 1)). Returns dict(difference, se, dof, t, low, high).

Parameters

Name	Type	Description	Default
mean_a	float	First sample mean, in a common measurement unit.	required
sd_a	float	Nonnegative sample standard deviation, in the mean measurement unit.	required
n_a	float	First sample count, at least two.	required
mean_b	float	Second sample mean, in the same unit as mean_a.	required
sd_b	float	Nonnegative sample standard deviation, in the mean measurement unit.	required
n_b	float	Second sample count, at least two.	required
confidence	float	Central confidence probability in (0, 1).	required

Returns

Name	Type	Description
	dict	Difference, se, low, high in input units; dof and t dimensionless.

Notes

Equation: SE = sqrt(sd_a^2/n_a+sd_b2/n_b); limits = mean_a-mean_b +/- t*SE.

Assumptions: Independent samples, at least two observations each, nonzero combined spread.

Excludes: Paired designs and exact guarantees for arbitrary heavy-tailed data.

References

Derivation: Welch (1947): use Satterthwaite effective degrees of freedom for unequal sample variances.

Examples

>>> calc_welch_interval(1, 1, 10, 1, 1, 10, 0.95)["difference"] == 0
True

calc_wilson_clean_runs

physics.statistics.calc_wilson_clean_runs(target_upper, z)

Find the clean-run count for a Wilson upper bound.

Fewest consecutive clean runs whose Wilson upper limit on the failure rate falls below target_upper.

With zero failures in n runs the upper limit is z^2 / (n + z^2), which is below t once n > z^2 (1 - t) / t.

Parameters

Name	Type	Description	Default
target_upper	float	Failure probability threshold in (0, 1), dimensionless.	required
z	float	Positive dimensionless normal quantile.	required

Returns

Name	Type	Description
	int	Run count making the upper bound less than the target.

Notes

Equation: n = floor(z^2*(1-target_upper)/target_upper)+1.

Assumptions: Independent Bernoulli failures; no failures observed; positive normal quantile.

Excludes: Optional-stopping guarantees and correlated failures.

References

Derivation: Solve the zero-success Wilson upper limit z^2/(n+z2) < target_upper.

Examples

>>> calc_wilson_clean_runs(0.1, 2) > calc_wilson_clean_runs(0.2, 2)
True

calc_wilson_interval

physics.statistics.calc_wilson_interval(successes, trials, z_critical)

Compute a Wilson score interval for a binomial proportion.

Parameters

Name	Type	Description	Default
successes	int	Accepted outcomes x, from zero through trials.	required
trials	int	Positive number n of independent Bernoulli trials.	required
z_critical	float	Nonnegative normal critical value; callers choose it for the desired coverage.	required

Returns

Name	Type	Description
	tuple of float	Lower and upper probability endpoints, clipped to [0, 1].

Notes

Equation: center=(x+z**2/2)/(n+z**2); half=z*sqrt(x*(n-x)/n+z**2/4)/(n+z**2).

Assumptions: Independent identically distributed Bernoulli trials and score-test inversion with a nonnegative critical value.

Excludes: Exact binomial coverage, clustering, and uncertainty from adaptive sample selection.

References

Wilson (1927), Probable Inference, the Law of Succession, and Statistical Inference. Derivation: Solve (p_hat-p)**2 <= z**2*p*(1-p)/n as a quadratic inequality in p.

Examples

>>> low, high = calc_wilson_interval(0, 10, 1.96)
>>> 0 <= low < high <= 1
True

calc_within_task_variance

physics.statistics.calc_within_task_variance(mean_rate, between_task_variance)

Recover within-task Bernoulli variance from task heterogeneity.

Mean within-task run variance E[p_i (1 - p_i)] = mean (1 - mean) - Var(p_i).

Parameters

Name	Type	Description	Default
mean_rate	float	Mean task probability in [0, 1], dimensionless.	required
between_task_variance	float	Nonnegative variance of task probabilities, dimensionless.	required

Returns

Name	Type	Description
	float	Dimensionless within-task outcome variance.

Notes

Equation: E[p_i(1-p_i)] = mean_rate*(1-mean_rate)-between_task_variance.

Assumptions: A valid mixture of Bernoulli task probabilities.

Excludes: Estimation uncertainty or correction of inconsistent variance inputs.

References

Derivation: Conditional on task i, a Bernoulli outcome has variance p_i*(1-p_i). Averaging over tasks gives E[p_i]-E[p_i^2]. Substitute E[p_i^2] = Var(p_i)+E[p_i]^2 to obtain the implemented decomposition.

Examples

>>> calc_within_task_variance(0.5, 0.0) == 0.25
True

fit_exponential_growth

physics.statistics.fit_exponential_growth(xs, ys, period=1.0)

Fit a constant multiplicative growth trend in log space.

Parameters

Name	Type	Description	Default
xs	sequence of float	At least two predictor coordinates in a common unit, with nonzero spread.	required
ys	sequence of float	Positive observations, equal in length to xs and expressed in a common unit.	required
period	float	Interval in x units for reporting growth; default 1.0.	1.0

Returns

Name	Type	Description
	tuple of float	Multiplicative factor per period, log10 slope per x unit, and log10 intercept.

Notes

Equation: Fit log10(y)=intercept+slope*x by ordinary least squares; factor=10**(slope*period).

Assumptions: Positive responses and a linear trend in logarithmic response; all points receive equal weight.

Excludes: Uncertainty intervals, heteroskedastic weights, structural breaks, and causal interpretation.

Returns (factor, slope, intercept) where log10(y) = intercept + slope * x and factor = 10 ** (slope * period). Trend lines on log-scale technology plots (compute, memory bandwidth, interconnect against year) are this fit.

References

Derivation: least-squares normal equations yield slope=Cov(x,log10(y))/Var(x), then exponentiation converts additive log growth to a multiplier.

Examples

>>> fit_exponential_growth([0.0, 1.0], [1.0, 10.0])
(10.0, 1.0, 0.0)