Univariate Statistics#

process_improve.univariate.metrics.t_value(p, v)[source]#

Return the value on the x-axis if you plot the cumulative t-distribution with a fractional area of p (p is therefore a fractional value between 0 and 1 on the y-axis) and v is the degrees of freedom.

Examples

Since the cumulative distribution passes symmetrically through the x-axis at 0.0 for any number of degrees of freedom

>>> t_value(0.5, v=10)
0.0

Zero fractional area under the curve is always at \(-\infty\):

>>> t_value(0.0, v=10)
-inf

100% fractional area is always at \(+\infty\):

>>> t_value(1.0, v=10)
inf

See also

t_value_cdf

does the inverse of this function.

Parameters:
Return type:

float

process_improve.univariate.metrics.t_value_cdf(z, v)[source]#

Return the fractional area under the cumulative t-distribution (y-axis value) at the t-value z on the x-axis, with v degrees of freedom.

Examples

The cumulative distribution is symmetric through the x-axis at 0.0 for any number of degrees of freedom, so half of the area lies below zero:

>>> t_value_cdf(0.0, v=10)
0.5

Zero fractional area under the curve is at \(-\infty\):

>>> t_value_cdf(float("-inf"), v=10)
0.0

100% fractional area is at \(+\infty\):

>>> t_value_cdf(float("inf"), v=10)
1.0

See also

t_value

does the inverse of this function.

Parameters:
Return type:

float

process_improve.univariate.metrics.test_normality(x)[source]#

Check the p-value of the hypothesis that the data are from a normal distribution.

If the p-value is less than the chosen alpha level (e.g. 0.05 or 0.025), then there is evidence that the data tested are NOT normally distributed.

On the other hand, if the p-value is greater than the chosen alpha level, then the null hypothesis that the data came from a normally distributed population can not be rejected. NOTE: it does not mean that the data are normally distributed, just that we have nothing better to say about it. See the Shapiro-Wilk test.

Implementation: Uses the Shapiro Wilk test directly taken from scipy.stats.shapiro.

Parameters:

x (ndarray | Series)

Return type:

float

process_improve.univariate.metrics.Sn(x, constant=1.1926)[source]#

Compute a robust scale estimator. The Sn metric is an efficient alternative to MAD.

Parameters:

  • x (np.ndarray or pd.Series) – A vector of values. NaN entries are ignored.

  • constant (float, optional) – Multiplicative constant that makes the estimator consistent with iid values from a Gaussian distribution with no outliers. Default is 1.1926.

Returns:

A scalar value, the Sn estimate of spread. Returns NaN if all entries are missing, and 0.0 when only a single non-missing value is supplied.

Return type:

np.floating

Notes

Follows the Rousseeuw-Croux definition [2] exactly, including the high median inside and the low median outside, and reproduces robustbase::Sn [1] (with finite.corr=TRUE) for both odd and even sample sizes.

Disadvantages of MAD:

  • It does not have particularly high efficiency for data that is in fact normal (37%). In comparison, the median has 64% efficiency for normal data.

  • The MAD statistic also has an implicit assumption of symmetry. That is, it measures the distance from a measure of central location (the median).

References

process_improve.univariate.metrics.ttest_independent(sample_A, sample_B, conflevel=0.995, *, equal_var=True)[source]#

Core calculation for a test of differences between the average of A and the average of B. No checking of inputs.

Parameters:

  • sample_A (iterable) – Vector of n_a measurements.

  • sample_B (iterable) – Vector of n_b measurements.

  • conflevel (float) – Value between 0 and 1 (closer to 1.0), that gives the level of confidence required for the 2-sided test.

  • equal_var (bool) – True (the default) runs Student’s pooled-variance test: the two samples are assumed to share one variance, which is estimated from both and spent on n_A + n_B - 2 degrees of freedom. False runs Welch’s test: each sample keeps its own variance and the degrees of freedom come from the Welch-Satterthwaite approximation. Matches scipy.stats.ttest_ind()’s parameter of the same name.

Returns:

Outcomes from the statistical test. The keys of interest are "t value" (the test statistic; its "p value" is computed from the t distribution on "Degrees of freedom"), "Std error of difference" (the denominator of that statistic, and the half-width unit of the confidence interval), and "Pooled std dev". "Equal variance assumed" records which test was run; "Pooled std dev" is NaN under Welch, where no pooled estimate exists.

Changed in version 2.0.0: The aliases "z value" and "Pooled standard deviation", deprecated since 1.70.0, are removed. Use "t value" and "Std error of difference".

Return type:

dict

Notes

The default stays Student’s so existing results do not move, but it is the weaker choice. Pooling is only valid when the two variances really are equal; when they are not, and especially when the larger variance sits with the smaller sample, Student’s test does not hold its nominal error rate, while Welch’s does and costs almost nothing when the variances are equal. Delacre, Lakens & Leys (2017) make the case that Welch’s should be the default, and it already is in R’s t.test: https://rips-irsp.com/articles/10.5334/irsp.82

Testing the variances first and choosing on the outcome is worse than either, not better: the pre-test spends its own error rate, and the conditional procedure has no clean size. Pick equal_var from what is known about the measurement, not from the data in hand.

process_improve.univariate.metrics.ttest_independent_from_df(df, grouper_column, values_column, conflevel=0.995, correction=None, *, equal_var=True)[source]#

Calculate the t-test for differences between two or more groups and returns a confidence interval for the difference. The test is for UNPAIRED differences.

The dataframe df contains a grouper_column with 2 or more unique values (e.g. ‘A’ and ‘B’). All unique values of the grouper_column are used, and t-tests are done between the values in the values_column.

Args:

df (pd.DataFrame): Dataframe of the values and grouping variable. grouper_column (str): Indicates which column will be grouped on. values_column (str): Which column contains the numeric values to calculate the test on. conflevel (float, optional): [description]. Defaults to 0.995. correction (str | None, optional): “holm”, “bh”, or None. Defaults to None. See below. equal_var (bool, optional): Forwarded to ttest_independent(). Defaults to True.

Equal variances: equal_var=True (the default) is Student’s pooled-variance test; False is Welch’s. Welch is the safer choice for a pairwise family, where the groups have no particular reason to share one variance; see ttest_independent()’s Notes.

Multiplicity: by default the returned p-values are UNCORRECTED. With k groups this runs k(k-1)/2 tests, so the chance of at least one false positive is well above the per-test level. Pass correction="holm" to control the family-wise error rate (see holm_bonferroni()) or correction="bh" to control the false discovery rate (see benjamini_hochberg()). Either adds p value (adjusted), reject and correction columns; the raw p value column is always kept.

Output: Dataframe with columns containing the statistical outputs of the t-test, including:

  1. Group “A” name

  2. Group “B” name

  3. Group “A” mean

  4. Group “B” mean

  5. t value for the difference between group “B” minus group “A”

  6. p-value for this t value

  7. Confidence interval low value for difference between group “B” minus group “A”

  8. Confidence interval high value for difference between group “B” minus group “A”

Example: df : has 3 levels in the grouper variable; [‘Marco’, ‘Pete’, ‘Sam’]

Output will have 3 rows: Group A name Group B name Marco Pete Marco Sam Pete Sam

Parameters:
Return type:

DataFrame

process_improve.univariate.metrics.ttest_paired(differences, conflevel=0.995)[source]#

Core calculation for a test of paired differences.

Parameters:

  • differences (pd.Series) – The paired differences (e.g. sample_A - sample_B) for which the test is run.

  • conflevel (float, optional) – Value between 0 and 1 (closer to 1.0), that gives the level of confidence required for the 2-sided test. Default is 0.995.

Returns:

Outcomes from the statistical test, with keys:

  • "Differences mean": the mean of the input differences.

  • "t value": the t statistic, the mean divided by its standard error.

  • "ConfInt: Lo", "ConfInt: Hi": lower and upper bounds of the two-sided confidence interval around the mean difference at the requested conflevel.

  • "p value": two-sided p-value from the t-distribution.

  • "Degrees of freedom": n - 1 for n paired observations.

  • "Std error of mean difference": s / sqrt(n), the scale the confidence interval is built on.

  • "Std dev of differences": s, the sample standard deviation of differences.

Deprecated since 2.1.0, to be removed in 3.0: "z value", the t statistic, and "Standard deviation", the standard error of the mean difference, under names that say otherwise. Reading either by name gives a DeprecationWarning.

Return type:

dict

process_improve.univariate.metrics.ttest_paired_from_df(df, grouper_column, values_column, conflevel=0.995, correction=None)[source]#

Calculate the t-test for paired differences between two or more groups and returns a confidence interval for the difference. The test is for PAIRED differences. The differences is always defined as the A values minus the B values: after - before, or A - B.

The dataframe df contains a grouper_column with 2 or more unique values (e.g. ‘A’ and ‘B’). All unique values of the grouper_column are used, and t-tests are done between the values in the values_column.

When selecting the columns, the number of values per column must be the same.

Args:

df (pd.DataFrame): Dataframe of the values and grouping variable. grouper_column (str): Indicates which column will be grouped on. values_column (str): Which column contains the numeric values to calculate the test on. conflevel (float, optional): [description]. Defaults to 0.995. correction (str | None, optional): “holm”, “bh”, or None. Defaults to None. See below.

Multiplicity: by default the returned p-values are UNCORRECTED. With k groups this runs k(k-1)/2 tests, so the chance of at least one false positive is well above the per-test level. Pass correction="holm" to control the family-wise error rate (see holm_bonferroni()) or correction="bh" to control the false discovery rate (see benjamini_hochberg()). Either adds p value (adjusted), reject and correction columns; the raw p value column is always kept.

Output: Dataframe with columns containing the statistical outputs of the t-test, including:

  1. Group A name

  2. Group B name

  3. Group A mean

  4. Group B mean

  5. Differences mean: average difference between the groups (not the same as the difference of the averages from items 3 and 4 above)

  6. t value for the difference between group “A” minus group “B”

  7. p value for this t value

  8. Confidence interval low value for difference between group “A” minus group “B”

  9. Confidence interval high value for difference between group “A” minus group “B”

  10. Std error of mean difference, and Std dev of differences

The columns z value and Standard deviation hold the t value and the standard error of the mean difference. They are deprecated since 2.1.0 and will be removed in 3.0; a column read cannot warn, so this note is the notice.

Parameters:
Return type:

DataFrame

process_improve.univariate.metrics.confidence_interval(df, column_name, conflevel=0.95, style='robust')[source]#

Calculate the confidence interval, returned as a tuple, for the column_name (str) in the dataframe df, for a given confidence level conflevel (default: 0.95).

style: [‘robust’; ‘regular’]: indicates which style of estimates to use for the center and

spread. Default: ‘robust’

Missing values are ignored.

Raises:

ValueError – If style is neither ‘robust’ nor ‘regular’, or if the column holds fewer than 2 non-missing values. An unrecognised style used to fall through to the classical branch, so a typo returned a different interval with nothing to signal it (#561).

Parameters:
Return type:

tuple

process_improve.univariate.metrics.median_absolute_deviation(x, axis=0, center=<function median>, scale='normal', nan_policy='omit')[source]#

Taken from scipy.stats.stats: we want the same functionality, but with a slightly different default function signature.

  • scale=’normal’ instead of scale=1.0.

  • nan_policy=’omit’ instead of nan_policy=’propogate’

Compute the median absolute deviation of the data along the given axis. The median absolute deviation (MAD) computes the median over the absolute deviations from the median. It is a measure of dispersion similar to the standard deviation but more robust to outliers. The MAD of an empty array is np.nan.

Parameters:

  • x (array_like) – Input array or object that can be converted to an array.

  • axis (int or None, optional) – Axis along which the range is computed. Default is 0. Axis == None will not be accepted anymore.

  • center (callable, optional) – A function that will return the central value. The default is to use np.median. Any user defined function used will need to have the function signature func(arr, axis).

  • scale (scalar or str, optional) – The numerical value of scale will be divided out of the final result. The default is "normal", which results in scale being the inverse of the standard normal quantile function at 0.75, which is approximately 0.67449 (so the returned value is consistent with the standard deviation for normally distributed data). Any positive float may also be passed. Array-like scale is also allowed, as long as it broadcasts correctly to the output such that out / scale is a valid operation. The output dimensions depend on the input array, x, and the axis argument.

  • nan_policy ({'propagate', 'raise', 'omit'}, optional) – Defines how to handle when input contains nan. The following options are available (default is ‘omit’): * ‘propagate’: returns nan * ‘raise’: throws an error * ‘omit’: performs the calculations ignoring nan values

Returns:

mad – If the input contains integers or floats of smaller precision than np.float64, then the output data-type is np.float64. Otherwise, the output data-type is the same as that of the input.

Return type:

scalar or ndarray

See also

numpy.std, numpy.var, numpy.median, scipy.stats.iqr, scipy.stats.tmean, scipy.stats.tstd, scipy.stats.tvar

Notes

The center argument only affects the calculation of the central value around which the MAD is calculated. That is, passing in center=np.mean will calculate the MAD around the mean - it will not calculate the mean absolute deviation.

The input array may contain inf, but if center returns inf, the corresponding MAD for that data will be nan.

References

Examples

When comparing the behavior of median_abs_deviation with np.std, the latter is affected when we change a single value of an array to have an outlier value while the MAD hardly changes:

>>> from scipy import stats
>>> x = stats.norm.rvs(size=100, scale=1, random_state=123456)
>>> x.std()
0.9973906394005013
>>> stats.median_abs_deviation(x)
0.82832610097857
>>> x[0] = 345.6
>>> x.std()
34.42304872314415
>>> stats.median_abs_deviation(x)
0.8323442311590675

Axis handling example:

>>> x = np.array([[10, 7, 4], [3, 2, 1]])
>>> x
array([[10,  7,  4],
       [ 3,  2,  1]])
>>> stats.median_abs_deviation(x)
array([3.5, 2.5, 1.5])
>>> stats.median_abs_deviation(x, axis=None)  <-- syntax of `axis=None` will be depricated
2.0

Scale normal example:

>>> x = stats.norm.rvs(size=1000000, scale=2, random_state=123456)
>>> stats.median_abs_deviation(x)
1.3487398527041636
>>> stats.median_abs_deviation(x, scale='normal')
1.9996446978061115
process_improve.univariate.metrics.biweight_midvariance(x, nan_policy='omit')[source]#

Return the Mosteller-Tukey robust scale (biweight midvariance) of x.

The biweight midvariance is a robust, highly efficient estimator of the variance: it down-weights observations far from the median and ignores gross outliers entirely.

Parameters:

  • x (np.ndarray or pd.Series) – One-dimensional sample of numeric values.

  • nan_policy ({"omit", "propagate"}, optional) – "omit" (default) drops missing values; "propagate" returns nan if any value is missing.

Returns:

The robust variance estimate. Returns 0.0 when the MAD is zero (e.g. a constant sample), and nan for an empty sample.

Return type:

float

References

Mosteller and Tukey, Data Analysis and Regression, pp. 207-208, 1977.

process_improve.univariate.metrics.holm_bonferroni(p_values, alpha=0.05)[source]#

Holm-Bonferroni step-down correction for multiple comparisons.

Holm’s method controls the family-wise error rate while being uniformly more powerful than the plain Bonferroni correction. It is the recommended post-hoc correction for a family of pairwise comparisons.

Parameters:

  • p_values (array-like) – The raw (uncorrected) p-values of the individual comparisons.

  • alpha (float, optional) – Family-wise significance level, by default 0.05.

Returns:

A bunch with, in the same order as the input:

  • p_adjusted: the Holm-adjusted p-values.

  • reject: boolean array, True where the null hypothesis is rejected at level alpha.

  • alpha: the family-wise level used.

Return type:

sklearn.utils.Bunch

References

Holm, “A simple sequentially rejective multiple test procedure”, Scandinavian Journal of Statistics, 6, 65-70, 1979.

process_improve.univariate.metrics.benjamini_hochberg(p_values, alpha=0.05)[source]#

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

The Benjamini-Hochberg (BH) procedure controls the expected proportion of false positives among the rejected hypotheses (the false discovery rate), rather than the family-wise error rate that holm_bonferroni() controls. BH is the conventional choice when a family contains many comparisons and a few false positives are tolerable, for example testing one product effect across many sensory attributes.

Parameters:

  • p_values (array-like) – The raw (uncorrected) p-values of the individual comparisons.

  • alpha (float, optional) – Target false-discovery rate, by default 0.05.

Returns:

A bunch with, in the same order as the input:

  • p_adjusted: the BH-adjusted p-values (q-values), monotone in the rank order of the input p-values.

  • reject: boolean array, True where the hypothesis is rejected at the target false-discovery rate alpha.

  • alpha: the target false-discovery rate used.

Return type:

sklearn.utils.Bunch

References

Benjamini and Hochberg, “Controlling the false discovery rate: a practical and powerful approach to multiple testing”, Journal of the Royal Statistical Society B, 57, 289-300, 1995.

process_improve.univariate.metrics.summary_stats(x, method='robust')[source]#

Return summary statistics of the numeric values in vector x.

Parameters:

  • x (numpy.ndarray or pandas.Series) – A vector of univariate values to summarize.

  • method (str, optional) – If "robust" (the default), the reported center is the median and the spread is the Sn robust estimate; otherwise the mean and the sample standard deviation are used.

Returns:

A summary of the univariate vector. The most useful keys are "center" (a measure of the center, e.g. the median for the robust method) and "spread" (a measure of the spread, e.g. the Sn robust estimate for the robust method).

Return type:

dict

process_improve.univariate.metrics.detect_outliers_esd(x, algorithm='esd', max_outliers_detected=1, **kwargs)[source]#

Return a list of indexes of points in the vector x which are likely outliers.

A second output (can be ignored) contains the details of the values used to make the decision.

Arguments:

x {list, sequence, NumPy vector/array} – [A sequence, list or vector which can be unravelled.]

Keyword Arguments:

algorithm – Two algorithms are possible to detect outliers: (default: “esd”)

‘esd’: Generalized ESD Test for Outliers. If max_outliers_detected=1 this is essentially Grubb’s test. For more details, please see: https://www.itl.nist.gov/div898/handbook/eda/section3/eda35h3.htm

‘cc-robust’: Build a robust control-chart for the sequence x and points should lie outside the +/- 3 sigma limits are considered outliers. Not Implemented Yet: left here as an idea for the future, but not confirmed yet. Passing algorithm='cc-robust' (or any other value besides 'esd') silently returns an empty outlier list and an empty details dict; it does not raise.

max_outliers_detected – The maximum number of outliers that

should be detected, as required by the algorithms.

kwargs – Algorithm dependent arguments. Defaults are shown here.

‘esd’:

‘robust_variant’ = False. When True, uses the median and MAD for the center and the spread respectively. WARNING: the ESD critical values (lambda_i) are derived for the classical statistic max|x - mean| / std, which is bounded above by (N-1)/sqrt(N); a MAD-scaled statistic has no such bound and is routinely 2-3x larger on contaminated data, so the robust variant declares far more outliers than the nominal alpha implies (it can flag points in perfectly clean data). It was previously the default; it is now opt-in and should be treated as a screening heuristic, not a calibrated test.

‘alpha’ = 0.05. The significance level of the testing.

Parameters:
Return type:

tuple[list[int], defaultdict[Any, Any]]

process_improve.univariate.metrics.tietjen_moore_test(x, n_outliers, *, two_sided=True, alpha=0.05, n_simulations=10000, random_state=None)[source]#

Tietjen-Moore test for a specified number of outliers.

Tests the null hypothesis “there are no outliers” against the alternative “the n_outliers most extreme observations are outliers”. Unlike the generalised ESD test, the number of suspected outliers must be fixed in advance. The test statistic has no closed-form critical value, so it is obtained by simulation under the normal null.

Parameters:

  • x (np.ndarray or pd.Series) – One-dimensional sample. Missing values are dropped.

  • n_outliers (int) – The number of suspected outliers to test for (1 <= n_outliers < N).

  • two_sided (bool, optional) – If True (default) the test looks for outliers on either tail (the observations with the largest absolute deviation from the mean). If False it tests only the n_outliers largest observations.

  • alpha (float, optional) – Significance level, by default 0.05.

  • n_simulations (int, optional) – Number of Monte-Carlo samples used to estimate the critical value.

  • random_state (int or None, optional) – Seed for the simulation, for reproducibility.

Returns:

With statistic, critical_value, reject (True when the outliers are significant), outlier_indices (positions in the missing-value-removed sample), n_outliers and alpha.

Return type:

sklearn.utils.Bunch

References

Tietjen and Moore, “Some Grubbs-type statistics for the detection of several outliers”, Technometrics, 14, 583-597, 1972. See also the NIST handbook, section 3.5.h.3.

process_improve.univariate.metrics.distribution_fit(x, distribution='norm', alpha=0.05)[source]#

Check how well a sample fits a named distribution.

Fits the parameters of the requested scipy.stats distribution by maximum likelihood and runs a Kolmogorov-Smirnov goodness-of-fit test (NIST handbook, section 3.5.7).

Parameters:

  • x (np.ndarray or pd.Series) – One-dimensional sample. Missing values are dropped.

  • distribution (str, optional) – Name of any continuous scipy.stats distribution, by default "norm".

  • alpha (float, optional) – Significance level for the fits_well verdict, by default 0.05.

Returns:

With distribution, fitted parameters, ks_statistic, ks_pvalue, fits_well (True when the fit is not rejected at level alpha) and the sample size n.

Return type:

sklearn.utils.Bunch

Notes

Because the distribution parameters are estimated from the same data, the KS p-value is conservative (the true Type-I error is smaller than alpha); it remains a useful screening check.

process_improve.univariate.metrics.variance_decomposition(df, measured, repeat)[source]#

Given a DataFrame df of raw data, and an indication of which column is the measured value column, and which is the repeat indicator, it will calculate the within and between replicate standard deviation.

Example Two measurements on day 1 [101, 102] and two measurements on day 2 [94, 95]. The between-day variation can already be expected to be much greater than the within-day variation.

>>> df = pd.DataFrame(data={'Result': [101, 102, 94, 95], 'Repeat': [1, 1, 2, 2]})
    Result  Repeat
0     101       1
1     102       1
2      94       2
3      95       2
>>> output = variance_decomposition(df, measured="Result", repeat="Repeat")
{'total_ms':       16.666667,
 'total_dof':      3,
 'within_ms':      0.5,
 'within_stddev':  0.70711,
 'within_dof':     2,
 'between_ms':     49.0,
 'between_stddev': 4.9244,
 'between_dof':    1}

Notes

  • SSQ = sum of squares

  • DOF = degrees of freedom

  • MS = mean square = (sum of squares) / (degrees of freedom) = SSQ / DOF = variance

  • between_ms is the raw ANOVA mean square, which estimates the sum of the within-group variance and n0 times the between-group variance (n0 = average group size), NOT the between-group variance itself. between_stddev reports the actual between-group variance COMPONENT, the square root of max(0, (MS_between - MS_within) / n0). Earlier versions reported sqrt(MS_between), which mixes the within-group variance into the “between” number.

Parameters:
Return type:

dict

Reference distributions from past records#

A reference distribution from a record of past results, to judge a new difference of averages.

A plant trial compares the average of \(n\) runs by one method with the average of the next \(n\) by another. Whether the difference is unusual can be judged against the same process’s own history, without assuming that the observations are independent or normally distributed (Box, Hunter and Hunter, chapter 2):

  • All differences. Every difference of two adjacent averages of \(n\) consecutive past observations (the later less the earlier) forms the reference set. The significance level of an observed difference is the share of the set at least as extreme. A member equal to it counts half, as Box, Hunter and Hunter count a tie.

  • Nearly independent differences. Adjacent members of that set share observations, so they are correlated. Taking only the differences of groups of \(2n\) a gap of one observation apart gives \(m\) nearly independent differences. Their root mean square about zero, \(s = \sqrt{\sum d^2/m}\), scales the observed difference as a t statistic on \(m\) degrees of freedom.

Examples

Box, Hunter and Hunter’s 210 past batch yields, against the plant trial’s difference of 1.30 between averages of 10:

>>> import numpy as np
>>> from process_improve.experiments.reference_catalogue import reference_experiments
>>> [record] = [e for e in reference_experiments(source="box-hunter-hunter-1978")
...             if e["name"] == "bhh-production-record-210"]
>>> history = [row[1] for row in record["runs"]["rows"]]
>>> found = reference_distribution(history, 10, observed=1.30)
>>> found["n_differences"], found["n_beyond"], round(found["significance"], 3)
(191, 9, 0.047)
process_improve.univariate.reference_distribution.ALTERNATIVES = ('greater', 'less')#

"greater" when the later average is expected to be larger, "less" when it is expected to be smaller.

Type:

How the observed difference is judged

process_improve.univariate.reference_distribution.TIE_TOLERANCE = 1e-09#

Differences closer than this share of the record’s spread to the observed one are ties.

process_improve.univariate.reference_distribution.reference_distribution(history, window, observed=None, alternative='greater', gap=1)[source]#

Build the reference distribution of differences of adjacent averages, and judge a difference by it.

Parameters:

  • history (Sequence[float]) – Past observations of the process, in time order.

  • window (int) – \(n\), the number of observations in each average: as many as each method had in the trial.

  • observed (float or None) – The trial’s difference, the later average less the earlier; None gives the reference distribution alone.

  • alternative (str) – "greater" (the later method is expected to give larger values) or "less".

  • gap (int) – The number of past observations skipped between the groups of \(2n\) that give the nearly independent differences.

Returns:

  • differences: every difference of adjacent averages, in time order, and n_differences, their number;

  • nearly_independent: the differences of groups of \(2n\) a gap apart, s their root mean square about zero, and df their number;

  • with observed: n_beyond and n_tied, the members beyond it and equal to it, significance (n_beyond plus half n_tied, over n_differences), t (the observed difference over s) and t_significance, the share of a t distribution on df degrees of freedom beyond t.

Return type:

dict[str, Any]

Raises:

ValueError – If the history is not finite numbers, window or gap is not a whole number of at least 1 (0 or more for gap), the history is too short for one difference, or alternative is not one of ALTERNATIVES.

Error transmission#

Error transmission: how the errors of measured quantities pass into a quantity worked out from them.

A quantity \(y = f(x_1, \dots, x_k)\) worked out from measured quantities carries their errors. When each input varies a little about its mean \(\mu_i\), with standard deviation \(\sigma_i\), and the errors are independent, the straight line that touches \(f\) at the means (its local linearization) gives the variance of \(y\):

\[\operatorname{Var}(y) \approx \sum_i \left(\frac{\partial f}{\partial x_i}\right)^2 \sigma_i^2 .\]

Each term is the variance that one input transmits, so the terms say which measurement to improve first (Box, Hunter and Hunter, section 17.2). The approximation holds while \(f\) is nearly straight over a few standard deviations of each input. The linearity check works \(f\) out at \(\mu_i - 3\sigma_i\) and \(\mu_i + 3\sigma_i\), the other inputs at their means, and sets each change beside the straight line’s \(\pm 3\sigma_i\,\partial f/\partial x_i\).

Examples

The volume of a sphere from its diameter, \(V = \pi d^3 / 6\), measured with a standard deviation of 0.01 about 2.0:

>>> import math
>>> from process_improve.univariate import error_transmission
>>> found = error_transmission(lambda d: math.pi * d**3 / 6, means={"d": 2.0}, std_devs={"d": 0.01})
>>> round(found["value"], 4), round(found["slopes"]["d"], 4), round(found["std_dev"], 4)
(4.1888, 6.2832, 0.0628)
process_improve.univariate.error_transmission.CHECK_SDS = 3.0#

How far each input moves for the linearity check, in its own standard deviations.

process_improve.univariate.error_transmission.STEP = 0.0001#

The step of the central difference that gives each slope, as a share of the input’s standard deviation.

process_improve.univariate.error_transmission.error_transmission(func, means, std_devs)[source]#

Return the variance a function transmits from the errors of its inputs, by local linearization.

Each slope \(\partial f/\partial x_i\) is a central difference at the means, over a step of STEP of the input’s standard deviation (or of its mean’s size, for an input known exactly). The transmitted variance is \(\sum_i (\partial f/\partial x_i)^2 \sigma_i^2\), for inputs whose errors are independent.

Parameters:

  • func (Callable[..., float]) – The quantity worked out from the inputs, called with each input by name: func(d=2.0, h=5.0).

  • means (Mapping[str, float]) – Each input’s mean: the point the function is linearized at.

  • std_devs (Mapping[str, float]) – Each input’s standard deviation, for the same inputs; 0 for an input known exactly.

Returns:

  • value: the function at the means.

  • variance and std_dev: the transmitted variance, and its square root.

  • slopes: each input’s partial derivative at the means.

  • contributions: the variance each input transmits, \((\partial f/\partial x_i)^2 \sigma_i^2\).

  • shares: each contribution as a fraction of the variance; empty when the variance is 0.

  • linearity: for each input with a standard deviation above 0, the function at \(\mu_i - 3\sigma_i\) and \(\mu_i + 3\sigma_i\) (low and high), the changes from value there (down and up), and the straight line’s change, \(3\sigma_i\,\partial f/\partial x_i\) (linear). When -down and up are both close to linear, the function is close to a straight line over the input’s spread, and the variance can be trusted.

Return type:

dict[str, Any]

Raises:

ValueError – If means and std_devs name different inputs, a mean is not a finite number, a standard deviation is negative or not a finite number, or the function is not a finite number at a point it is worked out at.

Examples

The variance of a product, \(y = x_1 x_2\), is close to \(x_2^2 \sigma_1^2 + x_1^2 \sigma_2^2\):

>>> found = error_transmission(lambda a, b: a * b, means={"a": 3.0, "b": 4.0}, std_devs={"a": 0.1, "b": 0.2})
>>> round(found["variance"], 6), {k: round(v, 6) for k, v in found["contributions"].items()}
(0.52, {'a': 0.16, 'b': 0.36})

Serial correlation#

Records whose neighbours in time are correlated: a first-order moving average, and the variance of a difference.

Successive batches of a process are often not independent: a high yield tends to follow a low one, or another high one. Two consequences matter when such a record is used as a reference.

  • The variance of a difference of averages. With a lag-1 autocorrelation \(\rho_1\) and none beyond, the variance of the difference of two adjacent averages of \(n\) consecutive observations is \((2\sigma^2/n)\,[1 + \rho_1 (2n - 3)/n]\), not the \(2\sigma^2/n\) of independent observations (Box, Hunter and Hunter, equation 18.5). A negative \(\rho_1\) shrinks it, so a difference judged as if the observations were independent looks less significant than it is.

  • A model of the record. The first-order moving average \(y_t = \mu + \varepsilon_t - \theta\varepsilon_{t-1}\) has a lag-1 autocorrelation \(-\theta/(1 + \theta^2)\) and none beyond it, which is what such a record shows.

Examples

Box, Hunter and Hunter’s 210 batch yields have \(s^2 = 8.30\) and \(r_1 = -0.29\); for two averages of 10 batches:

>>> from process_improve.univariate import difference_variance
>>> round(difference_variance(10, variance=8.30), 3), round(difference_variance(10, -0.29, 8.30), 3)
(1.66, 0.842)
process_improve.univariate.serial_correlation.MA1_METHODS = ('moments', 'likelihood')#

The ways fit_ma1() estimates the model.

process_improve.univariate.serial_correlation.difference_variance(n, rho1=0.0, variance=1.0, *, adjacent=True)[source]#

Return the variance of the difference of two averages of n consecutive observations.

The observations have variance \(\sigma^2\) and lag-1 autocorrelation \(\rho_1\), with none beyond lag 1 (a first-order moving average has that pattern). The variance of one average is \((\sigma^2/n)\,[1 + 2\rho_1(n - 1)/n]\). Two adjacent averages, the second block starting where the first ends, are themselves correlated through the one pair of neighbours that straddles them, so the variance of their difference is

\[\frac{2\sigma^2}{n}\left[1 + \rho_1 \frac{2n - 3}{n}\right],\]

which is Box, Hunter and Hunter’s \(1 + 1.7\rho_1\) for \(n = 10\) (equation 18.5). Two averages far apart are not correlated, and the bracket is \(1 + 2\rho_1(n - 1)/n\).

Parameters:

  • n (int) – The number of observations in each average.

  • rho1 (float) – The lag-1 autocorrelation of the observations, from -0.5 to 0.5: no process with autocorrelation at lag 1 only has one larger in size. 0 for independent observations.

  • variance (float) – The variance \(\sigma^2\) of one observation; 1 gives the variance in units of it.

  • adjacent (bool) – True when the second average’s observations follow the first’s directly, as in a reference set of successive averages; False for two averages far apart.

Returns:

The variance of the difference of the two averages.

Return type:

float

Raises:

ValueError – If n is not a whole number of at least 1, rho1 is outside -0.5 to 0.5, or the variance is negative or not a finite number.

process_improve.univariate.serial_correlation.fit_ma1(y, method='moments')[source]#

Fit a first-order moving average, \(y_t = \mu + \varepsilon_t - \theta\varepsilon_{t-1}\), to a record.

method="moments" matches the model to the record’s lag-1 autocorrelation \(r_1\) and variance \(s^2\), as Box, Hunter and Hunter do (section 18.2): the model’s lag-1 autocorrelation is \(-\theta/(1 + \theta^2)\), so \(\theta\) is the root of \(r_1\theta^2 + \theta + r_1 = 0\) between -1 and 1, and \(\sigma_\varepsilon^2 = s^2/(1 + \theta^2)\). No such root exists when \(r_1\) is 0.5 or more in size.

method="likelihood" maximizes the likelihood of the whole record (statsmodels’ ARIMA of order (0, 0, 1)), which uses every autocorrelation, not only the first, and gives a standard error for \(\theta\). The two estimates can differ a good deal on a short record.

Parameters:

  • y (Sequence[float]) – The record, in time order: at least three observations, not all the same.

  • method (str) – "moments" or "likelihood".

Returns:

mean (\(\hat\mu\)), theta, sigma2 (\(\hat\sigma_\varepsilon^2\)), r1 (the record’s lag-1 autocorrelation), theta_std_error (None for moments), n and method.

Return type:

dict[str, Any]

Raises:

ValueError – If the method is not one of MA1_METHODS, the record has fewer than three observations, a value is not a finite number, every value is the same, or (by moments) the lag-1 autocorrelation is 0.5 or more in size.

Examples

The first twelve of Box, Hunter and Hunter’s batch yields:

>>> found = fit_ma1([85.5, 81.7, 80.6, 84.7, 88.2, 84.9, 81.8, 84.9, 85.2, 81.9, 89.4, 79.0])
>>> round(found["r1"], 3), round(found["theta"], 3), round(found["sigma2"], 2)
(-0.337, 0.388, 8.19)

Results against a stated value#

Judge results against values stated before the data: a known mean, a new observation, a norm.

Box, Hunter and Hunter (chapters 2 and 6) set results against values that were known, or stated, before the data were taken:

  • A variance about a known mean. When the mean \(\eta\) is known, the squared deviations from it estimate \(\sigma^2\) on all \(n\) degrees of freedom, \(\sum (y - \eta)^2 / n\): no degree of freedom is spent estimating the mean.

  • One new observation. A single new value \(y\) is set against the stated mean and scaled by the standard deviation of earlier data, \(t = (y - \eta)/s\), on the degrees of freedom of \(s\): \(n\) when \(s\) is taken about the known mean, \(n - 1\) when about the earlier data’s own average.

  • A grand mean against a norm. In a one-way layout the sum of squares about a norm \(\eta_0\) splits into three parts,

    \[\sum_{ij} (y_{ij} - \eta_0)^2 = N(\bar y - \eta_0)^2 + \sum_i n_i (\bar y_i - \bar y)^2 + \sum_{ij} (y_{ij} - \bar y_i)^2,\]

    on 1, \(a - 1\) and \(N - a\) degrees of freedom, and the grand mean is tested against the norm with the residual’s mean square, which the differences between the groups do not inflate.

Examples

Box, Hunter and Hunter’s five cars (Exercise 2.12), whose mean is known to be 25.0, and a new car of the same kind that gives 28.2:

>>> cars = [27.8, 24.3, 22.8, 26.0, 24.2]
>>> round(variance_about(cars, 25.0)["variance"], 3)
2.962
>>> round(ttest_new_observation(28.2, 25.0, cars)["t_value"], 2)
1.86
process_improve.univariate.stated_values.variance_about(x, mean)[source]#

Return the variance of observations about a known mean, on all \(n\) degrees of freedom.

Parameters:

  • x (Sequence[float] or numpy.ndarray) – The observations.

  • mean (float) – The mean \(\eta\), known before the data: 0 for differences that average to zero, for example.

Returns:

variance (\(\sum (y - \eta)^2 / n\)), std_dev (its square root) and df (\(n\)).

Return type:

dict[str, float]

Raises:

ValueError – If there are no observations, or a value or the mean is not a finite number.

process_improve.univariate.stated_values.ttest_new_observation(y, mean, x, mean_known=True)[source]#

Test one new observation against a stated mean, scaled by the standard deviation of earlier data.

\(t = (y - \eta)/s\), with \(s\) from the earlier observations x: about the stated mean on \(n\) degrees of freedom when the mean is known, or about their own average on \(n - 1\) when it is not.

Parameters:

  • y (float) – The new observation.

  • mean (float) – The stated mean \(\eta\).

  • x (Sequence[float] or numpy.ndarray) – Earlier observations of the same kind, which give \(s\).

  • mean_known (bool) – Whether mean is the known mean of the earlier observations, so that \(s\) is taken about it.

Returns:

t_value, df, std_dev (\(s\)), the two-sided p_value and the one-sided p_value_greater and p_value_less.

Return type:

dict[str, float]

Raises:

ValueError – If a value is not a finite number, or x leaves no degree of freedom for \(s\).

process_improve.univariate.stated_values.grand_mean_test(y, groups, value)[source]#

Split a one-way layout’s sum of squares about a stated value, and test the grand mean against it.

Parameters:

  • y (Sequence[float] or numpy.ndarray) – The observations.

  • groups (Sequence[Hashable]) – The group of each observation (a treatment, a diet, a supplier).

  • value (float) – The norm \(\eta_0\), stated before the data.

Returns:

  • sum_of_squares and df, each with the keys total (about the norm), mean (\(N(\bar y - \eta_0)^2\), 1 degree of freedom), groups (between the group averages, \(a - 1\)) and residual (within the groups, \(N - a\)); the first is the sum of the other three;

  • grand_mean, \(\bar y\), and std_error, \(\sqrt{s^2 / N}\) with \(s^2\) the residual mean square;

  • t_value (\((\bar y - \eta_0)\) over its standard error) on df_resid (\(N - a\)) degrees of freedom, with the two-sided p_value and the one-sided p_value_greater and p_value_less.

Return type:

dict[str, Any]

Raises:

ValueError – If y and groups differ in length, a value is not a finite number, or the groups leave no residual degree of freedom (every group a single observation).

Examples

Box, Hunter and Hunter’s blood coagulation times of 24 animals on four diets (Table 6.1), against the norm of 63:

>>> times = [62, 60, 63, 59, 63, 67, 71, 64, 65, 66, 68, 66, 71, 67, 68, 68, 56, 62, 60, 61, 63, 64, 63, 59]
>>> diets = list("AAAA" "BBBBBB" "CCCCCC" "DDDDDDDD")
>>> found = grand_mean_test(times, diets, 63)
>>> {part: round(ss) for part, ss in found["sum_of_squares"].items()}
{'total': 364, 'mean': 24, 'groups': 228, 'residual': 112}
>>> round(found["t_value"], 2), found["df_resid"]
(2.07, 20)