mathematicskit.statistics#

mathematicskit.statistics hub · History · Examples · API reference

Descriptive statistics; hypothesis tests (z, t, chi-square, ANOVA); confidence intervals for means/proportions/variances; ordinary least-squares regression with residual diagnostics; and bootstrap resampling via scipy.stats.bootstrap.

mathematicskit.statistics: statistical inference.

Descriptive statistics (mean, variance, skewness, kurtosis, order statistics) via numpy/scipy.stats; hypothesis tests (one/two- sample z- and t-tests, chi-square goodness-of-fit and independence tests, one-way ANOVA) using mathematicskit.probability’s Normal for z-tests and scipy.stats’s t/chi-square/F distributions elsewhere; confidence intervals for means, proportions, and variances; ordinary least-squares linear regression with residual diagnostics and R^2/adjusted-R^2; bootstrap resampling via scipy.stats.bootstrap and the leave-one-out jackknife for confidence intervals, bias, and standard errors; Pearson and Spearman correlation; maximum-likelihood fitting and the likelihood-ratio test; distribution-free and exact tests (Kolmogorov-Smirnov, Fisher’s exact, Wilcoxon signed-rank, Mann-Whitney U); James-Stein shrinkage; and Bonferroni and Benjamini-Hochberg multiple-testing corrections.

class mathematicskit.statistics.BootstrapResult(estimate, lower, upper, std_error, confidence_level=0.95, method='BCa')[source]#

Bases: object

Container for a bootstrap confidence interval/standard error.

Parameters:
confidence_level: float = 0.95#

The interval’s nominal coverage probability.

Type:

float

estimate: float#

The statistic computed on the original (non-resampled) data.

Type:

float

lower: float#

Lower endpoint of the bootstrap confidence interval.

Type:

float

method: str = 'BCa'#

The bootstrap CI method (scipy.stats.bootstrap’s "percentile", "basic", or "BCa").

Type:

str

std_error: float#

The bootstrap distribution’s standard deviation – the resampling estimate of estimate’s standard error, obtained without any closed-form sampling-distribution assumption.

Type:

float

upper: float#

Upper endpoint of the bootstrap confidence interval.

Type:

float

class mathematicskit.statistics.ConfidenceIntervalResult(estimate, lower, upper, confidence_level=0.95, method='')[source]#

Bases: object

Container for a confidence interval.

Parameters:
confidence_level: float = 0.95#

The interval’s nominal coverage probability. The endpoints come from the \(\alpha/2\) and \(1-\alpha/2\) quantiles of the reference distribution, with \(\alpha = 1 - \text{confidence level}\).

Type:

float

estimate: float#

The point estimate (e.g. sample mean/proportion/variance).

Type:

float

lower: float#

The interval’s lower endpoint.

Type:

float

method: str = ''#

e.g. "z", "t", "chi_square_variance", "wald_proportion".

Type:

str

upper: float#

The interval’s upper endpoint.

Type:

float

width()[source]#

float: upper - lower.

Return type:

float

class mathematicskit.statistics.CorrelationResult(coefficient, p_value, n, method='')[source]#

Bases: object

Container for a correlation coefficient and its significance test.

Parameters:
coefficient: float#

The correlation coefficient, in \([-1, 1]\).

Type:

float

method: str = ''#

"pearson" or "spearman".

Type:

str

n: int#

Number of paired observations.

Type:

int

p_value: float#

Two-sided p-value for \(H_0\): no association.

Type:

float

class mathematicskit.statistics.DescriptiveStatsResult(n, mean, variance, std, skewness, kurtosis, minimum, q1, median, q3, maximum)[source]#

Bases: object

Container for a dataset’s descriptive-statistics summary.

Parameters:
kurtosis: float#

Sample excess kurtosis (scipy.stats.kurtosis, normal distribution has kurtosis 0 under this convention).

Type:

float

maximum: float#

The largest observation. Together with minimum, q1, median and q3 this makes up the five-number summary.

Type:

float

mean: float#

Sample mean.

Type:

float

median: float#

The 50th percentile.

Type:

float

minimum: float#

The smallest observation.

Type:

float

n: int#

Sample size.

Type:

int

q1: float#

First quartile (25th percentile).

Type:

float

q3: float#

Third quartile (75th percentile).

Type:

float

skewness: float#

Sample skewness (Fisher-Pearson, scipy.stats.skew).

Type:

float

std: float#

Sample standard deviation (ddof=1).

Type:

float

variance: float#

Sample variance (ddof=1, the unbiased estimator).

Type:

float

class mathematicskit.statistics.HypothesisTestResult(statistic, p_value, df=None, method='', extra=<factory>)[source]#

Bases: object

Container for the output of a hypothesis test.

Parameters:
df: float | None = None#

Degrees of freedom, where applicable.

Type:

float, optional

extra: dict#

Free-form diagnostics slot (e.g. group means for ANOVA).

Type:

dict

method: str = ''#

e.g. "one_sample_t", "chi_square_independence".

Type:

str

p_value: float#

The p-value under the null hypothesis.

Type:

float

reject_null(alpha=0.05)[source]#

Whether the p-value falls below the significance level alpha.

Parameters:

alpha (float)

Return type:

bool

Returns:

bool

statistic: float#

The test statistic (e.g. z, t, \(\chi^2\), F).

Type:

float

class mathematicskit.statistics.JackknifeResult(estimate, bias, std_error, bias_corrected, replicates)[source]#

Bases: object

Container for a jackknife bias and standard-error estimate.

Parameters:
bias: float#

Quenouille’s jackknife bias estimate.

Type:

float

bias_corrected: float#

estimate - bias.

Type:

float

estimate: float#

The statistic computed on the full sample.

Type:

float

replicates: ndarray#

The leave-one-out values of the statistic.

Type:

ndarray, shape (n,)

std_error: float#

Tukey’s jackknife standard-error estimate.

Type:

float

class mathematicskit.statistics.MaximumLikelihoodResult(params, log_likelihood, n_params, distribution='')[source]#

Bases: object

Container for a maximum-likelihood distribution fit.

Parameters:
aic()[source]#

float: Akaike’s information criterion, \(2k - 2\ell(\hat\theta)\).

Return type:

float

distribution: str = ''#

The scipy.stats distribution name, e.g. "norm".

Type:

str

log_likelihood: float#

The maximized log-likelihood \(\ell(\hat\theta)\).

Type:

float

n_params: int#

Number of parameters actually estimated (fixed ones excluded).

Type:

int

params: tuple#

The fitted parameters, in scipy.stats order (shape parameters first, then loc and scale).

Type:

tuple

class mathematicskit.statistics.MultipleTestingResult(rejected, adjusted_p_values, alpha=0.05, method='')[source]#

Bases: object

Container for a multiple-testing correction.

Parameters:
adjusted_p_values: ndarray#

Adjusted p-values, comparable directly with alpha.

Type:

ndarray, shape (m,)

alpha: float = 0.05#

The error rate being controlled (family-wise error rate or false discovery rate, depending on method).

Type:

float

method: str = ''#

"bonferroni" or "benjamini_hochberg".

Type:

str

rejected: ndarray#

Which null hypotheses are rejected.

Type:

ndarray of bool, shape (m,)

class mathematicskit.statistics.RegressionResult(coefficients, standard_errors, t_statistics, p_values, fitted_values, residuals, r_squared, adjusted_r_squared)[source]#

Bases: object

Container for an ordinary-least-squares regression fit.

Parameters:
adjusted_r_squared: float#

r_squared penalized for the number of predictors, \(1 - (1-R^2)(n-1)/(n-p)\). Unlike r_squared, it does not rise automatically each time another predictor is added, so it is the fairer figure when comparing models of different size.

Type:

float

coefficients: ndarray#

Fitted coefficients (intercept first, if included).

Type:

ndarray, shape (p,)

fitted_values: ndarray#

\(X\hat\beta\), the model’s prediction at each observed predictor.

Type:

ndarray, shape (n,)

p_values: ndarray#

Per-coefficient two-sided p-value.

Type:

ndarray, shape (p,)

r_squared: float#

The fraction of the response’s variance the fit explains, \(1 - SS_{\text{res}}/SS_{\text{tot}}\); 1.0 for a perfect fit.

Type:

float

residuals: ndarray#

y - fitted_values.

Type:

ndarray, shape (n,)

standard_errors: ndarray#

Standard error of each coefficient.

Type:

ndarray, shape (p,)

t_statistics: ndarray#

Per-coefficient t-statistic for H0: beta_j = 0.

Type:

ndarray, shape (p,)

mathematicskit.statistics.benjamini_hochberg(p_values, alpha=0.05)[source]#

Benjamini-Hochberg step-up procedure controlling the false discovery rate.

Sort the p-values \(p_{(1)} \le \dots \le p_{(m)}\), find the largest \(k\) with \(p_{(k)} \le k\alpha/m\), and reject \(H_{(1)}, \dots, H_{(k)}\). The adjusted p-values are \(\min_{j \ge i} \min(1, m p_{(j)}/j)\). See Benjamini & Hochberg (1995), J. R. Statist. Soc. B 57(1), 289-300.

Parameters:

  • p_values (ndarray)

  • alpha (float) – Target false discovery rate.

Return type:

MultipleTestingResult

Returns:

MultipleTestingResult

Examples

>>> result = benjamini_hochberg([0.01, 0.04, 0.03, 0.5])
>>> [round(float(v), 4) for v in result.adjusted_p_values]
[0.04, 0.0533, 0.0533, 0.5]
>>> result.rejected
array([ True, False, False, False])
mathematicskit.statistics.bonferroni_correction(p_values, alpha=0.05)[source]#

Bonferroni correction: reject \(H_i\) when \(m p_i \le \alpha\).

Controls the family-wise error rate (the chance of any false rejection) at \(\alpha\), under any dependence between tests.

Parameters:
Return type:

MultipleTestingResult

Returns:

MultipleTestingResult

Examples

>>> bonferroni_correction([0.01, 0.02, 0.2]).adjusted_p_values
array([0.03, 0.06, 0.6 ])
mathematicskit.statistics.bootstrap_confidence_interval(data, statistic=<function mean>, n_resamples=9999, confidence_level=0.95, method='BCa', seed=0)[source]#

Bootstrap confidence interval and standard error for an arbitrary statistic.

Resamples data with replacement n_resamples times, computes statistic on each resample, and builds a confidence interval from the resulting bootstrap distribution – the standard nonparametric alternative to a closed-form interval (e.g. mean_confidence_interval()) when no convenient parametric form is known, or as a check on one that is. See Efron & Tibshirani, An Introduction to the Bootstrap, 1993, Ch. 12-14.

Parameters:
Return type:

BootstrapResult

Returns:

BootstrapResult

Examples

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> data = rng.normal(loc=5.0, scale=2.0, size=200)
>>> result = bootstrap_confidence_interval(data, statistic=np.mean, seed=0)
>>> result.lower < 5.0 < result.upper
True
mathematicskit.statistics.chi_square_goodness_of_fit(observed, expected)[source]#

Chi-square goodness-of-fit test: does observed match expected proportions/counts?

\(\chi^2 = \sum_i (O_i - E_i)^2/E_i\), referred to scipy.stats.chi2 with \(k-1\) degrees of freedom (k categories).

The single degree of freedom subtracted is the one used up by the constraint that the counts sum to \(n\); this function assumes expected is fully specified in advance. If the expected counts were themselves obtained by estimating \(m\) parameters from the same data (fitting a Poisson rate, say), the correct degrees of freedom are \(k-1-m\), and the p-value reported here will be too large – conservative about rejecting \(H_0\). See DeGroot & Schervish, Probability and Statistics, 4th ed., Sec. 10.3.

Parameters:

  • observed (ndarray) – Observed counts.

  • expected (ndarray) – Expected counts under \(H_0\) (need not be pre-normalized to sum to the same total as observed, but should be positive).

Return type:

HypothesisTestResult

Returns:

HypothesisTestResult

Examples

>>> import numpy as np
>>> # A fair six-sided die rolled 600 times: expect 100 of each face.
>>> observed = np.array([90.0, 105.0, 95.0, 100.0, 110.0, 100.0])
>>> result = chi_square_goodness_of_fit(observed, expected=np.full(6, 100.0))
>>> result.df
5.0
>>> result.reject_null(alpha=0.05)
False
mathematicskit.statistics.chi_square_independence(contingency_table)[source]#

Chi-square test of independence for a two-way contingency table.

Expected counts under independence are \(E_{ij} = (\text{row total}_i)(\text{col total}_j)/N\); \(\chi^2 = \sum_{ij} (O_{ij}-E_{ij})^2/E_{ij}\), referred to scipy.stats.chi2 with \((r-1)(c-1)\) degrees of freedom. See DeGroot & Schervish, Probability and Statistics, 4th ed., Sec. 10.4.

Parameters:

contingency_table (ndarray) – Observed counts.

Return type:

HypothesisTestResult

Returns:

HypothesisTestResult – extra["expected"] holds the expected-counts table.

Examples

>>> import numpy as np
>>> table = np.array([[30.0, 10.0], [20.0, 40.0]])
>>> result = chi_square_independence(table)
>>> result.df
1.0
>>> result.reject_null(alpha=0.01)
True
mathematicskit.statistics.cohens_d(data1, data2)[source]#

Cohen’s d: a standardized (unit-free) mean difference between two groups.

\(d = \dfrac{\bar x_1 - \bar x_2}{s_{\text{pooled}}}\), where \(s_{\text{pooled}} = \sqrt{\dfrac{(n_1-1)s_1^2 + (n_2-1)s_2^2}{n_1+n_2-2}}\). Conventional (Cohen, 1988) rough benchmarks: \(|d|\approx 0.2\) small, \(0.5\) medium, \(0.8\) large – useful alongside two_sample_t_test(), since a test can be statistically significant (small p-value) while the underlying effect is practically negligible, especially at large sample sizes. See Cohen, Statistical Power Analysis for the Behavioral Sciences, 2nd ed., 1988, Ch. 2.

Parameters:
Return type:

float

Returns:

float

Examples

>>> import numpy as np
>>> a = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
>>> b = a + 1.0  # every value shifted up by 1
>>> round(cohens_d(a, b), 4)
-0.6325
mathematicskit.statistics.descriptive_stats(data)[source]#

Compute a standard descriptive-statistics summary of a dataset.

Mean, variance/standard deviation (ddof=1, the unbiased sample estimator), Fisher-Pearson skewness, excess kurtosis, and the five-number summary (min, quartiles, max). See DeGroot & Schervish, Probability and Statistics, 4th ed., Sec. 9.1-9.2.

Parameters:

data (ndarray)

Return type:

DescriptiveStatsResult

Returns:

DescriptiveStatsResult

Examples

>>> result = descriptive_stats([2.0, 4.0, 4.0, 4.0, 5.0, 5.0, 7.0, 9.0])
>>> result.n
8
>>> result.mean
5.0
>>> round(result.std, 4)
2.1381
mathematicskit.statistics.fisher_exact_test(table, alternative='two-sided')[source]#

Fisher’s exact test of independence for a 2x2 contingency table.

Conditioning on both margins, the top-left count follows a hypergeometric distribution under independence, so the p-value is an exact sum of hypergeometric probabilities rather than a chi-square approximation. Wraps scipy.stats.fisher_exact(). See R. A. Fisher, The Design of Experiments (1935), Ch. 2.

Parameters:

  • table (ndarray) – Observed counts.

  • alternative (str)

Return type:

HypothesisTestResult

Returns:

HypothesisTestResult – statistic is the sample odds ratio \(ad/bc\).

Examples

>>> # The lady tasting tea: all 8 cups classified correctly.
>>> result = fisher_exact_test([[4, 0], [0, 4]], alternative="greater")
>>> round(result.p_value * 70, 10)  # exactly 1/70
1.0
mathematicskit.statistics.jackknife(data, statistic=<function mean>)[source]#

Jackknife bias and standard-error estimates for a statistic.

Recomputes statistic on each of the \(n\) samples that leave one observation out. For the sample mean the jackknife standard error is exactly \(s/\sqrt n\); for the plug-in variance \(\frac1n\sum(x_i-\bar x)^2\) the bias-corrected estimate is exactly the unbiased \(\frac1{n-1}\sum(x_i-\bar x)^2\).

Parameters:
Return type:

JackknifeResult

Returns:

JackknifeResult

Examples

>>> import numpy as np
>>> data = np.array([2.0, 4.0, 4.0, 5.0, 7.0, 9.0])
>>> result = jackknife(data, np.var)  # plug-in variance, ddof=0
>>> bool(np.isclose(result.bias_corrected, np.var(data, ddof=1)))
True
mathematicskit.statistics.james_stein_estimator(x, sigma=1.0, target=None, positive_part=True)[source]#

James-Stein estimate of a normal mean vector.

Shrinks x toward target (the origin by default) by the factor \(1 - (p-2)\sigma^2/\lVert x - \text{target}\rVert^2\). The positive-part version clips that factor at 0, which never does worse and avoids overshooting past the target.

Parameters:

  • x (ndarray) – One observation of each of \(p \ge 3\) normal means.

  • sigma (float) – Known common standard deviation of each observation.

  • target (ndarray | None) – Point to shrink toward; defaults to the origin.

  • positive_part (bool) – Clip the shrinkage factor at 0.

Return type:

ndarray

Returns:

ndarray, shape (p,)

Examples

>>> import numpy as np
>>> james_stein_estimator(np.array([2.0, 2.0, 2.0, 2.0]))  # factor 1 - 2/16
array([1.75, 1.75, 1.75, 1.75])
mathematicskit.statistics.kolmogorov_smirnov_test(data, reference='norm', args=(), alternative='two-sided')[source]#

Kolmogorov-Smirnov test against a distribution or a second sample.

The statistic is the largest vertical gap between the empirical CDF \(F_n\) and the reference CDF \(F\), \(D_n = \sup_x |F_n(x) - F(x)|\). Kolmogorov (1933) showed that \(\sqrt n D_n\) has the same limiting distribution for every continuous \(F\); Smirnov (1939) extended it to two samples. Wraps scipy.stats.kstest() (one sample) or scipy.stats.ks_2samp() (two samples).

Parameters:

  • data (ndarray)

  • reference (str | Callable | ndarray) – A scipy.stats distribution name or a CDF callable (one-sample test), or a second sample (two-sample test).

  • args (tuple) – Parameters of the reference distribution, e.g. (loc, scale). They must be fixed in advance, not estimated from data.

  • alternative (str)

Return type:

HypothesisTestResult

Returns:

HypothesisTestResult

Examples

>>> import numpy as np
>>> result = kolmogorov_smirnov_test(np.array([0.1, 0.4, 0.7]), "uniform")
>>> round(result.statistic, 6)  # largest gap is at x = 0.7: 1 - 0.7
0.3
mathematicskit.statistics.likelihood_ratio_test(log_likelihood_null, log_likelihood_alt, df)[source]#

Likelihood-ratio test of a nested null model against a larger model.

The statistic \(\Lambda = 2(\ell_1 - \ell_0)\) is, by Wilks’s theorem, asymptotically \(\chi^2_{df}\) under \(H_0\), where df is the number of extra free parameters in the alternative model. See Wilks (1938), Ann. Math. Statist. 9(1), 60-62.

Parameters:

  • log_likelihood_null (float) – Maximized log-likelihood of the restricted (null) model.

  • log_likelihood_alt (float) – Maximized log-likelihood of the full (alternative) model.

  • df (int) – Difference in the number of free parameters.

Return type:

HypothesisTestResult

Returns:

HypothesisTestResult

Examples

>>> result = likelihood_ratio_test(-10.0, -8.0, df=1)
>>> result.statistic
4.0
>>> round(result.p_value, 4)
0.0455
mathematicskit.statistics.linear_regression(x, y, add_intercept=True)[source]#

Ordinary least squares: \(y = X\beta + \varepsilon\).

Solves for \(\hat\beta\) via numpy.linalg.lstsq(); standard errors come from the diagonal of \(\hat\sigma^2 (X^TX)^{-1}\) where \(\hat\sigma^2 = SSR/(n-p)\), and each coefficient’s t-statistic/p-value tests \(H_0: \beta_j = 0\) against scipy.stats.t with \(n-p\) degrees of freedom. See DeGroot & Schervish, Probability and Statistics, 4th ed., Sec. 11.5-11.7.

Parameters:

  • x (ndarray) – Predictor(s). A 1D array is treated as a single predictor.

  • y (ndarray) – Response.

  • add_intercept (bool) – Whether to prepend a column of ones to x (fitting an intercept term). coefficients[0] is then the intercept.

Return type:

RegressionResult

Returns:

RegressionResult

Examples

>>> import numpy as np
>>> x = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
>>> y = 2.0 * x + 1.0  # noiseless: fit should be essentially exact
>>> result = linear_regression(x, y)
>>> np.allclose(result.coefficients, [1.0, 2.0], atol=1e-8)
True
>>> round(result.r_squared, 6)
1.0
mathematicskit.statistics.mann_whitney_u_test(x, y, alternative='two-sided')[source]#

Mann-Whitney U (Wilcoxon rank-sum) test for two independent samples.

\(U\) counts the pairs \((x_i, y_j)\) with \(x_i > y_j\) (ties count one half); \(U/(n_1 n_2)\) estimates \(P(X > Y)\). Wraps scipy.stats.mannwhitneyu(). See Wilcoxon (1945) and H. B. Mann & D. R. Whitney, Ann. Math. Statist. 18(1) (1947), 50-60.

Parameters:
Return type:

HypothesisTestResult

Returns:

HypothesisTestResult – statistic is \(U\) for x.

Examples

>>> result = mann_whitney_u_test([5.0, 6.0, 7.0], [1.0, 2.0, 3.0])
>>> result.statistic  # every x beats every y: U = 3 * 3
9.0
mathematicskit.statistics.maximum_likelihood_fit(data, distribution='norm', **fixed)[source]#

Fit a continuous distribution to data by maximum likelihood.

Finds \(\hat\theta = \arg\max_\theta \sum_i \log f(x_i;\theta)\) using scipy.stats.<distribution>.fit. For the normal distribution the maximizer is closed-form: \(\hat\mu = \bar x\) and \(\hat\sigma^2 = \frac1n\sum_i (x_i-\bar x)^2\) (note the \(1/n\), not \(1/(n-1)\)). See Fisher (1922), Phil. Trans. R. Soc. A 222, 309-368.

Parameters:

  • data (ndarray)

  • distribution (str | rv_continuous) – A continuous scipy.stats distribution or its name ("norm", "expon", "gamma", …).

  • **fixed – Parameters held fixed rather than estimated, using scipy’s f-prefix convention (e.g. floc=0).

Return type:

MaximumLikelihoodResult

Returns:

MaximumLikelihoodResult

Examples

>>> import numpy as np
>>> data = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
>>> result = maximum_likelihood_fit(data, "norm")
>>> [round(float(v), 6) for v in result.params]  # mean and sqrt(2)
[3.0, 1.414214]
mathematicskit.statistics.mean_confidence_interval(data, sigma=None, confidence_level=0.95)[source]#

Confidence interval for a population mean.

\(\bar x \pm z_{\alpha/2}\,\sigma/\sqrt n\) if sigma (the population standard deviation) is known, else \(\bar x \pm t_{\alpha/2,\,n-1}\,s/\sqrt n\) using the sample standard deviation. See DeGroot & Schervish, Probability and Statistics, 4th ed., Sec. 8.5-8.6.

Parameters:

  • data (ndarray)

  • sigma (float | None) – Known population standard deviation; if omitted, the sample standard deviation and Student’s t distribution are used instead.

  • confidence_level (float)

Return type:

ConfidenceIntervalResult

Returns:

ConfidenceIntervalResult

Examples

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> data = rng.normal(loc=10.0, scale=2.0, size=100)
>>> result = mean_confidence_interval(data, sigma=2.0)
>>> bool(result.lower < 10.0 < result.upper)
True
mathematicskit.statistics.one_sample_t_test(data, mu0, alternative='two-sided')[source]#

One-sample t-test for the mean, with unknown population standard deviation.

\(t = \dfrac{\bar x - \mu_0}{s/\sqrt n}\), referred to Student’s t distribution with \(n-1\) degrees of freedom (scipy.stats.t). See DeGroot & Schervish, Probability and Statistics, 4th ed., Sec. 9.5.

Parameters:
Return type:

HypothesisTestResult

Returns:

HypothesisTestResult

Examples

>>> import numpy as np
>>> rng = np.random.default_rng(1)
>>> data = rng.normal(loc=100.0, scale=15.0, size=40)
>>> result = one_sample_t_test(data, mu0=100.0)
>>> result.df
39.0
mathematicskit.statistics.one_sample_z_test(data, mu0, sigma, alternative='two-sided')[source]#

One-sample z-test for the mean, with known population standard deviation.

\(z = \dfrac{\bar x - \mu_0}{\sigma/\sqrt n}\), referred to mathematicskit.probability.Normal(0, 1). See DeGroot & Schervish, Probability and Statistics, 4th ed., Sec. 9.5.

Parameters:

  • data (ndarray)

  • mu0 (float) – Hypothesized mean under \(H_0\).

  • sigma (float) – Known population standard deviation.

  • alternative (str)

Return type:

HypothesisTestResult

Returns:

HypothesisTestResult

Examples

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> data = rng.normal(loc=6.0, scale=2.0, size=500)
>>> result = one_sample_z_test(data, mu0=5.0, sigma=2.0)
>>> bool(result.reject_null(alpha=0.05))
True
mathematicskit.statistics.one_way_anova(*groups)[source]#

One-way analysis of variance: do 2+ groups share a common mean?

\(F = \dfrac{\text{between-group variance}}{\text{within-group variance}} = \dfrac{SSB/(k-1)}{SSW/(N-k)}\), referred to scipy.stats.f with \((k-1, N-k)\) degrees of freedom. See DeGroot & Schervish, Probability and Statistics, 4th ed., Sec. 11.3.

Parameters:

*groups (ndarray) – Two or more samples (one per group).

Return type:

HypothesisTestResult

Returns:

HypothesisTestResult – extra["group_means"] holds each group’s sample mean.

Examples

>>> import numpy as np
>>> a = np.array([4.0, 5.0, 6.0, 5.0])
>>> b = np.array([7.0, 8.0, 7.0, 9.0])
>>> c = np.array([4.0, 3.0, 5.0, 4.0])
>>> result = one_way_anova(a, b, c)
>>> result.reject_null(alpha=0.01)
True
mathematicskit.statistics.order_statistic(data, k)[source]#

The k-th order statistic (1-indexed: k=1 is the minimum).

Parameters:
Return type:

float

Returns:

float

Examples

>>> order_statistic([5.0, 1.0, 3.0, 2.0, 4.0], k=1)
1.0
>>> order_statistic([5.0, 1.0, 3.0, 2.0, 4.0], k=5)
5.0
mathematicskit.statistics.pearson_correlation(x, y)[source]#

Pearson’s product-moment correlation coefficient.

\[r = \frac{\sum_i (x_i - \bar x)(y_i - \bar y)} {\sqrt{\sum_i (x_i - \bar x)^2 \sum_i (y_i - \bar y)^2}},\]

with the p-value from \(t = r\sqrt{(n-2)/(1-r^2)}\) referred to Student’s t with \(n-2\) degrees of freedom (exact for bivariate normal data). Wraps scipy.stats.pearsonr(). See Pearson (1896), Phil. Trans. R. Soc. A 187, 253-318.

Parameters:

  • x (ndarray) – Paired observations, \(n \ge 3\).

  • y (ndarray) – Paired observations, \(n \ge 3\).

Return type:

CorrelationResult

Returns:

CorrelationResult

Examples

>>> import numpy as np
>>> x = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
>>> result = pearson_correlation(x, 2.0 * x + 1.0)
>>> round(result.coefficient, 12)
1.0
mathematicskit.statistics.proportion_confidence_interval(successes, n, confidence_level=0.95)[source]#

Wald (normal-approximation) confidence interval for a population proportion.

\(\hat p \pm z_{\alpha/2}\sqrt{\hat p(1-\hat p)/n}\), valid for n large enough that the normal approximation to the binomial holds (a common rule of thumb: \(n\hat p \geq 5\) and \(n(1-\hat p) \geq 5\)). See DeGroot & Schervish, Probability and Statistics, 4th ed., Sec. 9.5.

Parameters:
Return type:

ConfidenceIntervalResult

Returns:

ConfidenceIntervalResult

Examples

>>> result = proportion_confidence_interval(successes=520, n=1000)
>>> round(result.estimate, 4)
0.52
mathematicskit.statistics.spearman_correlation(x, y)[source]#

Spearman’s rank correlation coefficient.

Pearson’s \(r\) computed on the ranks of x and y (ties get their average rank). Without ties it reduces to Spearman’s formula

\[\rho = 1 - \frac{6 \sum_i d_i^2}{n(n^2 - 1)},\]

where \(d_i\) is the difference between the two ranks of observation \(i\). Wraps scipy.stats.spearmanr(). See Spearman (1904), Am. J. Psychol. 15(1), 72-101.

Parameters:

  • x (ndarray) – Paired observations, \(n \ge 3\).

  • y (ndarray) – Paired observations, \(n \ge 3\).

Return type:

CorrelationResult

Returns:

CorrelationResult

Examples

>>> import numpy as np
>>> x = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
>>> round(spearman_correlation(x, np.exp(x)).coefficient, 12)  # any increasing map gives 1
1.0
mathematicskit.statistics.two_sample_t_test(data1, data2, equal_var=True, alternative='two-sided')[source]#

Two-sample t-test for a difference of means.

Pooled-variance (“Student’s”) t-test if equal_var=True, otherwise Welch’s t-test with the Welch-Satterthwaite degrees of freedom. Referred to scipy.stats.t. See DeGroot & Schervish, Probability and Statistics, 4th ed., Sec. 9.6, and Welch (1947), Biometrika 34.

Parameters:
Return type:

HypothesisTestResult

Returns:

HypothesisTestResult

Examples

>>> import numpy as np
>>> rng = np.random.default_rng(2)
>>> a = rng.normal(loc=10.0, scale=2.0, size=50)
>>> b = rng.normal(loc=11.0, scale=2.0, size=50)
>>> result = two_sample_t_test(a, b)
>>> bool(result.statistic < 0)
True
mathematicskit.statistics.two_sample_z_test(data1, data2, sigma1, sigma2, alternative='two-sided')[source]#

Two-sample z-test for a difference of means, with known population standard deviations.

\(z = \dfrac{\bar x_1 - \bar x_2}{\sqrt{\sigma_1^2/n_1 + \sigma_2^2/n_2}}\). See DeGroot & Schervish, Probability and Statistics, 4th ed., Sec. 9.6.

Parameters:

  • data1 (ndarray)

  • data2 (ndarray)

  • sigma1 (float) – Known population standard deviations.

  • sigma2 (float) – Known population standard deviations.

  • alternative (str)

Return type:

HypothesisTestResult

Returns:

HypothesisTestResult

mathematicskit.statistics.variance_confidence_interval(data, confidence_level=0.95)[source]#

Confidence interval for a population variance, assuming normality.

\(\left(\dfrac{(n-1)s^2}{\chi^2_{\alpha/2,\,n-1}}, \dfrac{(n-1)s^2}{\chi^2_{1-\alpha/2,\,n-1}}\right)\), via scipy.stats.chi2. See DeGroot & Schervish, Probability and Statistics, 4th ed., Sec. 9.8.

Parameters:
Return type:

ConfidenceIntervalResult

Returns:

ConfidenceIntervalResult

Examples

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> data = rng.normal(loc=0.0, scale=3.0, size=200)
>>> result = variance_confidence_interval(data)
>>> bool(result.lower < 9.0 < result.upper)
True
mathematicskit.statistics.wilcoxon_signed_rank_test(x, y=None, alternative='two-sided')[source]#

Wilcoxon signed-rank test for a median of zero (or paired differences).

Ranks the absolute differences \(|x_i - y_i|\) and sums the ranks of the positive ones, \(W^+\). Under \(H_0\) (differences symmetric about 0) each sign is a fair coin flip, which gives an exact null distribution. Wraps scipy.stats.wilcoxon(). See F. Wilcoxon, Biometrics Bulletin 1(6) (1945), 80-83.

Parameters:

  • x (ndarray)

  • y (ndarray | None) – Paired sample; if given, the test is applied to x - y.

  • alternative (str)

Return type:

HypothesisTestResult

Returns:

HypothesisTestResult – statistic is scipy’s convention: \(\min(W^+, W^-)\) for two-sided tests, \(W^+\) otherwise.

Examples

>>> # 5 positive differences: exact one-sided p = (1/2)^5
>>> result = wilcoxon_signed_rank_test([1.0, 2.0, 3.0, 4.0, 5.0], alternative="greater")
>>> result.p_value
0.03125
On this page
Show Source