finalizare 1.0
This commit is contained in:
andrei-mihnea-cerbu
@@ -0,0 +1,643 @@
|
||||
"""
|
||||
.. _statsrefmanual:
|
||||
|
||||
==========================================
|
||||
Statistical functions (:mod:`scipy.stats`)
|
||||
==========================================
|
||||
|
||||
.. currentmodule:: scipy.stats
|
||||
|
||||
This module contains a large number of probability distributions,
|
||||
summary and frequency statistics, correlation functions and statistical
|
||||
tests, masked statistics, kernel density estimation, quasi-Monte Carlo
|
||||
functionality, and more.
|
||||
|
||||
Statistics is a very large area, and there are topics that are out of scope
|
||||
for SciPy and are covered by other packages. Some of the most important ones
|
||||
are:
|
||||
|
||||
- `statsmodels <https://www.statsmodels.org/stable/index.html>`__:
|
||||
regression, linear models, time series analysis, extensions to topics
|
||||
also covered by ``scipy.stats``.
|
||||
- `Pandas <https://pandas.pydata.org/>`__: tabular data, time series
|
||||
functionality, interfaces to other statistical languages.
|
||||
- `PyMC <https://docs.pymc.io/>`__: Bayesian statistical
|
||||
modeling, probabilistic machine learning.
|
||||
- `scikit-learn <https://scikit-learn.org/>`__: classification, regression,
|
||||
model selection.
|
||||
- `Seaborn <https://seaborn.pydata.org/>`__: statistical data visualization.
|
||||
- `rpy2 <https://rpy2.github.io/>`__: Python to R bridge.
|
||||
|
||||
|
||||
Probability distributions
|
||||
=========================
|
||||
|
||||
Each univariate distribution is an instance of a subclass of `rv_continuous`
|
||||
(`rv_discrete` for discrete distributions):
|
||||
|
||||
.. autosummary::
|
||||
:toctree: generated/
|
||||
|
||||
rv_continuous
|
||||
rv_discrete
|
||||
rv_histogram
|
||||
|
||||
Continuous distributions
|
||||
------------------------
|
||||
|
||||
.. autosummary::
|
||||
:toctree: generated/
|
||||
|
||||
alpha -- Alpha
|
||||
anglit -- Anglit
|
||||
arcsine -- Arcsine
|
||||
argus -- Argus
|
||||
beta -- Beta
|
||||
betaprime -- Beta Prime
|
||||
bradford -- Bradford
|
||||
burr -- Burr (Type III)
|
||||
burr12 -- Burr (Type XII)
|
||||
cauchy -- Cauchy
|
||||
chi -- Chi
|
||||
chi2 -- Chi-squared
|
||||
cosine -- Cosine
|
||||
crystalball -- Crystalball
|
||||
dgamma -- Double Gamma
|
||||
dweibull -- Double Weibull
|
||||
erlang -- Erlang
|
||||
expon -- Exponential
|
||||
exponnorm -- Exponentially Modified Normal
|
||||
exponweib -- Exponentiated Weibull
|
||||
exponpow -- Exponential Power
|
||||
f -- F (Snecdor F)
|
||||
fatiguelife -- Fatigue Life (Birnbaum-Saunders)
|
||||
fisk -- Fisk
|
||||
foldcauchy -- Folded Cauchy
|
||||
foldnorm -- Folded Normal
|
||||
genlogistic -- Generalized Logistic
|
||||
gennorm -- Generalized normal
|
||||
genpareto -- Generalized Pareto
|
||||
genexpon -- Generalized Exponential
|
||||
genextreme -- Generalized Extreme Value
|
||||
gausshyper -- Gauss Hypergeometric
|
||||
gamma -- Gamma
|
||||
gengamma -- Generalized gamma
|
||||
genhalflogistic -- Generalized Half Logistic
|
||||
genhyperbolic -- Generalized Hyperbolic
|
||||
geninvgauss -- Generalized Inverse Gaussian
|
||||
gibrat -- Gibrat
|
||||
gompertz -- Gompertz (Truncated Gumbel)
|
||||
gumbel_r -- Right Sided Gumbel, Log-Weibull, Fisher-Tippett, Extreme Value Type I
|
||||
gumbel_l -- Left Sided Gumbel, etc.
|
||||
halfcauchy -- Half Cauchy
|
||||
halflogistic -- Half Logistic
|
||||
halfnorm -- Half Normal
|
||||
halfgennorm -- Generalized Half Normal
|
||||
hypsecant -- Hyperbolic Secant
|
||||
invgamma -- Inverse Gamma
|
||||
invgauss -- Inverse Gaussian
|
||||
invweibull -- Inverse Weibull
|
||||
jf_skew_t -- Jones and Faddy Skew-T
|
||||
johnsonsb -- Johnson SB
|
||||
johnsonsu -- Johnson SU
|
||||
kappa4 -- Kappa 4 parameter
|
||||
kappa3 -- Kappa 3 parameter
|
||||
ksone -- Distribution of Kolmogorov-Smirnov one-sided test statistic
|
||||
kstwo -- Distribution of Kolmogorov-Smirnov two-sided test statistic
|
||||
kstwobign -- Limiting Distribution of scaled Kolmogorov-Smirnov two-sided test statistic.
|
||||
laplace -- Laplace
|
||||
laplace_asymmetric -- Asymmetric Laplace
|
||||
levy -- Levy
|
||||
levy_l
|
||||
levy_stable
|
||||
logistic -- Logistic
|
||||
loggamma -- Log-Gamma
|
||||
loglaplace -- Log-Laplace (Log Double Exponential)
|
||||
lognorm -- Log-Normal
|
||||
loguniform -- Log-Uniform
|
||||
lomax -- Lomax (Pareto of the second kind)
|
||||
maxwell -- Maxwell
|
||||
mielke -- Mielke's Beta-Kappa
|
||||
moyal -- Moyal
|
||||
nakagami -- Nakagami
|
||||
ncx2 -- Non-central chi-squared
|
||||
ncf -- Non-central F
|
||||
nct -- Non-central Student's T
|
||||
norm -- Normal (Gaussian)
|
||||
norminvgauss -- Normal Inverse Gaussian
|
||||
pareto -- Pareto
|
||||
pearson3 -- Pearson type III
|
||||
powerlaw -- Power-function
|
||||
powerlognorm -- Power log normal
|
||||
powernorm -- Power normal
|
||||
rdist -- R-distribution
|
||||
rayleigh -- Rayleigh
|
||||
rel_breitwigner -- Relativistic Breit-Wigner
|
||||
rice -- Rice
|
||||
recipinvgauss -- Reciprocal Inverse Gaussian
|
||||
semicircular -- Semicircular
|
||||
skewcauchy -- Skew Cauchy
|
||||
skewnorm -- Skew normal
|
||||
studentized_range -- Studentized Range
|
||||
t -- Student's T
|
||||
trapezoid -- Trapezoidal
|
||||
triang -- Triangular
|
||||
truncexpon -- Truncated Exponential
|
||||
truncnorm -- Truncated Normal
|
||||
truncpareto -- Truncated Pareto
|
||||
truncweibull_min -- Truncated minimum Weibull distribution
|
||||
tukeylambda -- Tukey-Lambda
|
||||
uniform -- Uniform
|
||||
vonmises -- Von-Mises (Circular)
|
||||
vonmises_line -- Von-Mises (Line)
|
||||
wald -- Wald
|
||||
weibull_min -- Minimum Weibull (see Frechet)
|
||||
weibull_max -- Maximum Weibull (see Frechet)
|
||||
wrapcauchy -- Wrapped Cauchy
|
||||
|
||||
The ``fit`` method of the univariate continuous distributions uses
|
||||
maximum likelihood estimation to fit the distribution to a data set.
|
||||
The ``fit`` method can accept regular data or *censored data*.
|
||||
Censored data is represented with instances of the `CensoredData`
|
||||
class.
|
||||
|
||||
.. autosummary::
|
||||
:toctree: generated/
|
||||
|
||||
CensoredData
|
||||
|
||||
|
||||
Multivariate distributions
|
||||
--------------------------
|
||||
|
||||
.. autosummary::
|
||||
:toctree: generated/
|
||||
|
||||
multivariate_normal -- Multivariate normal distribution
|
||||
matrix_normal -- Matrix normal distribution
|
||||
dirichlet -- Dirichlet
|
||||
dirichlet_multinomial -- Dirichlet multinomial distribution
|
||||
wishart -- Wishart
|
||||
invwishart -- Inverse Wishart
|
||||
multinomial -- Multinomial distribution
|
||||
special_ortho_group -- SO(N) group
|
||||
ortho_group -- O(N) group
|
||||
unitary_group -- U(N) group
|
||||
random_correlation -- random correlation matrices
|
||||
multivariate_t -- Multivariate t-distribution
|
||||
multivariate_hypergeom -- Multivariate hypergeometric distribution
|
||||
random_table -- Distribution of random tables with given marginals
|
||||
uniform_direction -- Uniform distribution on S(N-1)
|
||||
vonmises_fisher -- Von Mises-Fisher distribution
|
||||
|
||||
`scipy.stats.multivariate_normal` methods accept instances
|
||||
of the following class to represent the covariance.
|
||||
|
||||
.. autosummary::
|
||||
:toctree: generated/
|
||||
|
||||
Covariance -- Representation of a covariance matrix
|
||||
|
||||
|
||||
Discrete distributions
|
||||
----------------------
|
||||
|
||||
.. autosummary::
|
||||
:toctree: generated/
|
||||
|
||||
bernoulli -- Bernoulli
|
||||
betabinom -- Beta-Binomial
|
||||
betanbinom -- Beta-Negative Binomial
|
||||
binom -- Binomial
|
||||
boltzmann -- Boltzmann (Truncated Discrete Exponential)
|
||||
dlaplace -- Discrete Laplacian
|
||||
geom -- Geometric
|
||||
hypergeom -- Hypergeometric
|
||||
logser -- Logarithmic (Log-Series, Series)
|
||||
nbinom -- Negative Binomial
|
||||
nchypergeom_fisher -- Fisher's Noncentral Hypergeometric
|
||||
nchypergeom_wallenius -- Wallenius's Noncentral Hypergeometric
|
||||
nhypergeom -- Negative Hypergeometric
|
||||
planck -- Planck (Discrete Exponential)
|
||||
poisson -- Poisson
|
||||
randint -- Discrete Uniform
|
||||
skellam -- Skellam
|
||||
yulesimon -- Yule-Simon
|
||||
zipf -- Zipf (Zeta)
|
||||
zipfian -- Zipfian
|
||||
|
||||
|
||||
An overview of statistical functions is given below. Many of these functions
|
||||
have a similar version in `scipy.stats.mstats` which work for masked arrays.
|
||||
|
||||
Summary statistics
|
||||
==================
|
||||
|
||||
.. autosummary::
|
||||
:toctree: generated/
|
||||
|
||||
describe -- Descriptive statistics
|
||||
gmean -- Geometric mean
|
||||
hmean -- Harmonic mean
|
||||
pmean -- Power mean
|
||||
kurtosis -- Fisher or Pearson kurtosis
|
||||
mode -- Modal value
|
||||
moment -- Central moment
|
||||
expectile -- Expectile
|
||||
skew -- Skewness
|
||||
kstat --
|
||||
kstatvar --
|
||||
tmean -- Truncated arithmetic mean
|
||||
tvar -- Truncated variance
|
||||
tmin --
|
||||
tmax --
|
||||
tstd --
|
||||
tsem --
|
||||
variation -- Coefficient of variation
|
||||
find_repeats
|
||||
rankdata
|
||||
tiecorrect
|
||||
trim_mean
|
||||
gstd -- Geometric Standard Deviation
|
||||
iqr
|
||||
sem
|
||||
bayes_mvs
|
||||
mvsdist
|
||||
entropy
|
||||
differential_entropy
|
||||
median_abs_deviation
|
||||
|
||||
Frequency statistics
|
||||
====================
|
||||
|
||||
.. autosummary::
|
||||
:toctree: generated/
|
||||
|
||||
cumfreq
|
||||
percentileofscore
|
||||
scoreatpercentile
|
||||
relfreq
|
||||
|
||||
.. autosummary::
|
||||
:toctree: generated/
|
||||
|
||||
binned_statistic -- Compute a binned statistic for a set of data.
|
||||
binned_statistic_2d -- Compute a 2-D binned statistic for a set of data.
|
||||
binned_statistic_dd -- Compute a d-D binned statistic for a set of data.
|
||||
|
||||
Hypothesis Tests and related functions
|
||||
======================================
|
||||
SciPy has many functions for performing hypothesis tests that return a
|
||||
test statistic and a p-value, and several of them return confidence intervals
|
||||
and/or other related information.
|
||||
|
||||
The headings below are based on common uses of the functions within, but due to
|
||||
the wide variety of statistical procedures, any attempt at coarse-grained
|
||||
categorization will be imperfect. Also, note that tests within the same heading
|
||||
are not interchangeable in general (e.g. many have different distributional
|
||||
assumptions).
|
||||
|
||||
One Sample Tests / Paired Sample Tests
|
||||
--------------------------------------
|
||||
One sample tests are typically used to assess whether a single sample was
|
||||
drawn from a specified distribution or a distribution with specified properties
|
||||
(e.g. zero mean).
|
||||
|
||||
.. autosummary::
|
||||
:toctree: generated/
|
||||
|
||||
ttest_1samp
|
||||
binomtest
|
||||
quantile_test
|
||||
skewtest
|
||||
kurtosistest
|
||||
normaltest
|
||||
jarque_bera
|
||||
shapiro
|
||||
anderson
|
||||
cramervonmises
|
||||
ks_1samp
|
||||
goodness_of_fit
|
||||
chisquare
|
||||
power_divergence
|
||||
|
||||
Paired sample tests are often used to assess whether two samples were drawn
|
||||
from the same distribution; they differ from the independent sample tests below
|
||||
in that each observation in one sample is treated as paired with a
|
||||
closely-related observation in the other sample (e.g. when environmental
|
||||
factors are controlled between observations within a pair but not among pairs).
|
||||
They can also be interpreted or used as one-sample tests (e.g. tests on the
|
||||
mean or median of *differences* between paired observations).
|
||||
|
||||
.. autosummary::
|
||||
:toctree: generated/
|
||||
|
||||
ttest_rel
|
||||
wilcoxon
|
||||
|
||||
Association/Correlation Tests
|
||||
-----------------------------
|
||||
|
||||
These tests are often used to assess whether there is a relationship (e.g.
|
||||
linear) between paired observations in multiple samples or among the
|
||||
coordinates of multivariate observations.
|
||||
|
||||
.. autosummary::
|
||||
:toctree: generated/
|
||||
|
||||
linregress
|
||||
pearsonr
|
||||
spearmanr
|
||||
pointbiserialr
|
||||
kendalltau
|
||||
weightedtau
|
||||
somersd
|
||||
siegelslopes
|
||||
theilslopes
|
||||
page_trend_test
|
||||
multiscale_graphcorr
|
||||
|
||||
These association tests and are to work with samples in the form of contingency
|
||||
tables. Supporting functions are available in `scipy.stats.contingency`.
|
||||
|
||||
.. autosummary::
|
||||
:toctree: generated/
|
||||
|
||||
chi2_contingency
|
||||
fisher_exact
|
||||
barnard_exact
|
||||
boschloo_exact
|
||||
|
||||
Independent Sample Tests
|
||||
------------------------
|
||||
Independent sample tests are typically used to assess whether multiple samples
|
||||
were independently drawn from the same distribution or different distributions
|
||||
with a shared property (e.g. equal means).
|
||||
|
||||
Some tests are specifically for comparing two samples.
|
||||
|
||||
.. autosummary::
|
||||
:toctree: generated/
|
||||
|
||||
ttest_ind_from_stats
|
||||
poisson_means_test
|
||||
ttest_ind
|
||||
mannwhitneyu
|
||||
bws_test
|
||||
ranksums
|
||||
brunnermunzel
|
||||
mood
|
||||
ansari
|
||||
cramervonmises_2samp
|
||||
epps_singleton_2samp
|
||||
ks_2samp
|
||||
kstest
|
||||
|
||||
Others are generalized to multiple samples.
|
||||
|
||||
.. autosummary::
|
||||
:toctree: generated/
|
||||
|
||||
f_oneway
|
||||
tukey_hsd
|
||||
dunnett
|
||||
kruskal
|
||||
alexandergovern
|
||||
fligner
|
||||
levene
|
||||
bartlett
|
||||
median_test
|
||||
friedmanchisquare
|
||||
anderson_ksamp
|
||||
|
||||
Resampling and Monte Carlo Methods
|
||||
----------------------------------
|
||||
The following functions can reproduce the p-value and confidence interval
|
||||
results of most of the functions above, and often produce accurate results in a
|
||||
wider variety of conditions. They can also be used to perform hypothesis tests
|
||||
and generate confidence intervals for custom statistics. This flexibility comes
|
||||
at the cost of greater computational requirements and stochastic results.
|
||||
|
||||
.. autosummary::
|
||||
:toctree: generated/
|
||||
|
||||
monte_carlo_test
|
||||
permutation_test
|
||||
bootstrap
|
||||
|
||||
Instances of the following object can be passed into some hypothesis test
|
||||
functions to perform a resampling or Monte Carlo version of the hypothesis
|
||||
test.
|
||||
|
||||
.. autosummary::
|
||||
:toctree: generated/
|
||||
|
||||
MonteCarloMethod
|
||||
PermutationMethod
|
||||
BootstrapMethod
|
||||
|
||||
Multiple Hypothesis Testing and Meta-Analysis
|
||||
---------------------------------------------
|
||||
These functions are for assessing the results of individual tests as a whole.
|
||||
Functions for performing specific multiple hypothesis tests (e.g. post hoc
|
||||
tests) are listed above.
|
||||
|
||||
.. autosummary::
|
||||
:toctree: generated/
|
||||
|
||||
combine_pvalues
|
||||
false_discovery_control
|
||||
|
||||
|
||||
The following functions are related to the tests above but do not belong in the
|
||||
above categories.
|
||||
|
||||
Quasi-Monte Carlo
|
||||
=================
|
||||
|
||||
.. toctree::
|
||||
:maxdepth: 4
|
||||
|
||||
stats.qmc
|
||||
|
||||
Contingency Tables
|
||||
==================
|
||||
|
||||
.. toctree::
|
||||
:maxdepth: 4
|
||||
|
||||
stats.contingency
|
||||
|
||||
Masked statistics functions
|
||||
===========================
|
||||
|
||||
.. toctree::
|
||||
|
||||
stats.mstats
|
||||
|
||||
|
||||
Other statistical functionality
|
||||
===============================
|
||||
|
||||
Transformations
|
||||
---------------
|
||||
|
||||
.. autosummary::
|
||||
:toctree: generated/
|
||||
|
||||
boxcox
|
||||
boxcox_normmax
|
||||
boxcox_llf
|
||||
yeojohnson
|
||||
yeojohnson_normmax
|
||||
yeojohnson_llf
|
||||
obrientransform
|
||||
sigmaclip
|
||||
trimboth
|
||||
trim1
|
||||
zmap
|
||||
zscore
|
||||
gzscore
|
||||
|
||||
Statistical distances
|
||||
---------------------
|
||||
|
||||
.. autosummary::
|
||||
:toctree: generated/
|
||||
|
||||
wasserstein_distance
|
||||
wasserstein_distance_nd
|
||||
energy_distance
|
||||
|
||||
Sampling
|
||||
--------
|
||||
|
||||
.. toctree::
|
||||
:maxdepth: 4
|
||||
|
||||
stats.sampling
|
||||
|
||||
Random variate generation / CDF Inversion
|
||||
-----------------------------------------
|
||||
|
||||
.. autosummary::
|
||||
:toctree: generated/
|
||||
|
||||
rvs_ratio_uniforms
|
||||
|
||||
Fitting / Survival Analysis
|
||||
---------------------------
|
||||
|
||||
.. autosummary::
|
||||
:toctree: generated/
|
||||
|
||||
fit
|
||||
ecdf
|
||||
logrank
|
||||
|
||||
Directional statistical functions
|
||||
---------------------------------
|
||||
|
||||
.. autosummary::
|
||||
:toctree: generated/
|
||||
|
||||
directional_stats
|
||||
circmean
|
||||
circvar
|
||||
circstd
|
||||
|
||||
Sensitivity Analysis
|
||||
--------------------
|
||||
|
||||
.. autosummary::
|
||||
:toctree: generated/
|
||||
|
||||
sobol_indices
|
||||
|
||||
Plot-tests
|
||||
----------
|
||||
|
||||
.. autosummary::
|
||||
:toctree: generated/
|
||||
|
||||
ppcc_max
|
||||
ppcc_plot
|
||||
probplot
|
||||
boxcox_normplot
|
||||
yeojohnson_normplot
|
||||
|
||||
Univariate and multivariate kernel density estimation
|
||||
-----------------------------------------------------
|
||||
|
||||
.. autosummary::
|
||||
:toctree: generated/
|
||||
|
||||
gaussian_kde
|
||||
|
||||
Warnings / Errors used in :mod:`scipy.stats`
|
||||
--------------------------------------------
|
||||
|
||||
.. autosummary::
|
||||
:toctree: generated/
|
||||
|
||||
DegenerateDataWarning
|
||||
ConstantInputWarning
|
||||
NearConstantInputWarning
|
||||
FitError
|
||||
|
||||
Result classes used in :mod:`scipy.stats`
|
||||
-----------------------------------------
|
||||
|
||||
.. warning::
|
||||
|
||||
These classes are private, but they are included here because instances
|
||||
of them are returned by other statistical functions. User import and
|
||||
instantiation is not supported.
|
||||
|
||||
.. toctree::
|
||||
:maxdepth: 2
|
||||
|
||||
stats._result_classes
|
||||
|
||||
""" # noqa: E501
|
||||
|
||||
from ._warnings_errors import (ConstantInputWarning, NearConstantInputWarning,
|
||||
DegenerateDataWarning, FitError)
|
||||
from ._stats_py import *
|
||||
from ._variation import variation
|
||||
from .distributions import *
|
||||
from ._morestats import *
|
||||
from ._multicomp import *
|
||||
from ._binomtest import binomtest
|
||||
from ._binned_statistic import *
|
||||
from ._kde import gaussian_kde
|
||||
from . import mstats
|
||||
from . import qmc
|
||||
from ._multivariate import *
|
||||
from . import contingency
|
||||
from .contingency import chi2_contingency
|
||||
from ._censored_data import CensoredData
|
||||
from ._resampling import (bootstrap, monte_carlo_test, permutation_test,
|
||||
MonteCarloMethod, PermutationMethod, BootstrapMethod)
|
||||
from ._entropy import *
|
||||
from ._hypotests import *
|
||||
from ._rvs_sampling import rvs_ratio_uniforms
|
||||
from ._page_trend_test import page_trend_test
|
||||
from ._mannwhitneyu import mannwhitneyu
|
||||
from ._bws_test import bws_test
|
||||
from ._fit import fit, goodness_of_fit
|
||||
from ._covariance import Covariance
|
||||
from ._sensitivity_analysis import *
|
||||
from ._survival import *
|
||||
|
||||
# Deprecated namespaces, to be removed in v2.0.0
|
||||
from . import (
|
||||
biasedurn, kde, morestats, mstats_basic, mstats_extras, mvn, stats
|
||||
)
|
||||
|
||||
|
||||
__all__ = [s for s in dir() if not s.startswith("_")] # Remove dunders.
|
||||
|
||||
from scipy._lib._testutils import PytestTester
|
||||
test = PytestTester(__name__)
|
||||
del PytestTester
|
||||
Binary file not shown.
BIN
Binary file not shown.
BIN
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
BIN
Binary file not shown.
Binary file not shown.
Binary file not shown.
BIN
Binary file not shown.
BIN
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
BIN
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
BIN
Binary file not shown.
Binary file not shown.
Binary file not shown.
BIN
Binary file not shown.
BIN
Binary file not shown.
Binary file not shown.
Binary file not shown.
BIN
Binary file not shown.
Binary file not shown.
BIN
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
BIN
Binary file not shown.
BIN
Binary file not shown.
@@ -0,0 +1,642 @@
|
||||
# Many scipy.stats functions support `axis` and `nan_policy` parameters.
|
||||
# When the two are combined, it can be tricky to get all the behavior just
|
||||
# right. This file contains utility functions useful for scipy.stats functions
|
||||
# that support `axis` and `nan_policy`, including a decorator that
|
||||
# automatically adds `axis` and `nan_policy` arguments to a function.
|
||||
|
||||
import numpy as np
|
||||
from functools import wraps
|
||||
from scipy._lib._docscrape import FunctionDoc, Parameter
|
||||
from scipy._lib._util import _contains_nan, AxisError, _get_nan
|
||||
import inspect
|
||||
|
||||
|
||||
def _broadcast_arrays(arrays, axis=None):
|
||||
"""
|
||||
Broadcast shapes of arrays, ignoring incompatibility of specified axes
|
||||
"""
|
||||
new_shapes = _broadcast_array_shapes(arrays, axis=axis)
|
||||
if axis is None:
|
||||
new_shapes = [new_shapes]*len(arrays)
|
||||
return [np.broadcast_to(array, new_shape)
|
||||
for array, new_shape in zip(arrays, new_shapes)]
|
||||
|
||||
|
||||
def _broadcast_array_shapes(arrays, axis=None):
|
||||
"""
|
||||
Broadcast shapes of arrays, ignoring incompatibility of specified axes
|
||||
"""
|
||||
shapes = [np.asarray(arr).shape for arr in arrays]
|
||||
return _broadcast_shapes(shapes, axis)
|
||||
|
||||
|
||||
def _broadcast_shapes(shapes, axis=None):
|
||||
"""
|
||||
Broadcast shapes, ignoring incompatibility of specified axes
|
||||
"""
|
||||
if not shapes:
|
||||
return shapes
|
||||
|
||||
# input validation
|
||||
if axis is not None:
|
||||
axis = np.atleast_1d(axis)
|
||||
axis_int = axis.astype(int)
|
||||
if not np.array_equal(axis_int, axis):
|
||||
raise AxisError('`axis` must be an integer, a '
|
||||
'tuple of integers, or `None`.')
|
||||
axis = axis_int
|
||||
|
||||
# First, ensure all shapes have same number of dimensions by prepending 1s.
|
||||
n_dims = max([len(shape) for shape in shapes])
|
||||
new_shapes = np.ones((len(shapes), n_dims), dtype=int)
|
||||
for row, shape in zip(new_shapes, shapes):
|
||||
row[len(row)-len(shape):] = shape # can't use negative indices (-0:)
|
||||
|
||||
# Remove the shape elements of the axes to be ignored, but remember them.
|
||||
if axis is not None:
|
||||
axis[axis < 0] = n_dims + axis[axis < 0]
|
||||
axis = np.sort(axis)
|
||||
if axis[-1] >= n_dims or axis[0] < 0:
|
||||
message = (f"`axis` is out of bounds "
|
||||
f"for array of dimension {n_dims}")
|
||||
raise AxisError(message)
|
||||
|
||||
if len(np.unique(axis)) != len(axis):
|
||||
raise AxisError("`axis` must contain only distinct elements")
|
||||
|
||||
removed_shapes = new_shapes[:, axis]
|
||||
new_shapes = np.delete(new_shapes, axis, axis=1)
|
||||
|
||||
# If arrays are broadcastable, shape elements that are 1 may be replaced
|
||||
# with a corresponding non-1 shape element. Assuming arrays are
|
||||
# broadcastable, that final shape element can be found with:
|
||||
new_shape = np.max(new_shapes, axis=0)
|
||||
# except in case of an empty array:
|
||||
new_shape *= new_shapes.all(axis=0)
|
||||
|
||||
# Among all arrays, there can only be one unique non-1 shape element.
|
||||
# Therefore, if any non-1 shape element does not match what we found
|
||||
# above, the arrays must not be broadcastable after all.
|
||||
if np.any(~((new_shapes == 1) | (new_shapes == new_shape))):
|
||||
raise ValueError("Array shapes are incompatible for broadcasting.")
|
||||
|
||||
if axis is not None:
|
||||
# Add back the shape elements that were ignored
|
||||
new_axis = axis - np.arange(len(axis))
|
||||
new_shapes = [tuple(np.insert(new_shape, new_axis, removed_shape))
|
||||
for removed_shape in removed_shapes]
|
||||
return new_shapes
|
||||
else:
|
||||
return tuple(new_shape)
|
||||
|
||||
|
||||
def _broadcast_array_shapes_remove_axis(arrays, axis=None):
|
||||
"""
|
||||
Broadcast shapes of arrays, dropping specified axes
|
||||
|
||||
Given a sequence of arrays `arrays` and an integer or tuple `axis`, find
|
||||
the shape of the broadcast result after consuming/dropping `axis`.
|
||||
In other words, return output shape of a typical hypothesis test on
|
||||
`arrays` vectorized along `axis`.
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> import numpy as np
|
||||
>>> from scipy.stats._axis_nan_policy import _broadcast_array_shapes
|
||||
>>> a = np.zeros((5, 2, 1))
|
||||
>>> b = np.zeros((9, 3))
|
||||
>>> _broadcast_array_shapes((a, b), 1)
|
||||
(5, 3)
|
||||
"""
|
||||
# Note that here, `axis=None` means do not consume/drop any axes - _not_
|
||||
# ravel arrays before broadcasting.
|
||||
shapes = [arr.shape for arr in arrays]
|
||||
return _broadcast_shapes_remove_axis(shapes, axis)
|
||||
|
||||
|
||||
def _broadcast_shapes_remove_axis(shapes, axis=None):
|
||||
"""
|
||||
Broadcast shapes, dropping specified axes
|
||||
|
||||
Same as _broadcast_array_shapes, but given a sequence
|
||||
of array shapes `shapes` instead of the arrays themselves.
|
||||
"""
|
||||
shapes = _broadcast_shapes(shapes, axis)
|
||||
shape = shapes[0]
|
||||
if axis is not None:
|
||||
shape = np.delete(shape, axis)
|
||||
return tuple(shape)
|
||||
|
||||
|
||||
def _broadcast_concatenate(arrays, axis, paired=False):
|
||||
"""Concatenate arrays along an axis with broadcasting."""
|
||||
arrays = _broadcast_arrays(arrays, axis if not paired else None)
|
||||
res = np.concatenate(arrays, axis=axis)
|
||||
return res
|
||||
|
||||
|
||||
# TODO: add support for `axis` tuples
|
||||
def _remove_nans(samples, paired):
|
||||
"Remove nans from paired or unpaired 1D samples"
|
||||
# potential optimization: don't copy arrays that don't contain nans
|
||||
if not paired:
|
||||
return [sample[~np.isnan(sample)] for sample in samples]
|
||||
|
||||
# for paired samples, we need to remove the whole pair when any part
|
||||
# has a nan
|
||||
nans = np.isnan(samples[0])
|
||||
for sample in samples[1:]:
|
||||
nans = nans | np.isnan(sample)
|
||||
not_nans = ~nans
|
||||
return [sample[not_nans] for sample in samples]
|
||||
|
||||
|
||||
def _remove_sentinel(samples, paired, sentinel):
|
||||
"Remove sentinel values from paired or unpaired 1D samples"
|
||||
# could consolidate with `_remove_nans`, but it's not quite as simple as
|
||||
# passing `sentinel=np.nan` because `(np.nan == np.nan) is False`
|
||||
|
||||
# potential optimization: don't copy arrays that don't contain sentinel
|
||||
if not paired:
|
||||
return [sample[sample != sentinel] for sample in samples]
|
||||
|
||||
# for paired samples, we need to remove the whole pair when any part
|
||||
# has a nan
|
||||
sentinels = (samples[0] == sentinel)
|
||||
for sample in samples[1:]:
|
||||
sentinels = sentinels | (sample == sentinel)
|
||||
not_sentinels = ~sentinels
|
||||
return [sample[not_sentinels] for sample in samples]
|
||||
|
||||
|
||||
def _masked_arrays_2_sentinel_arrays(samples):
|
||||
# masked arrays in `samples` are converted to regular arrays, and values
|
||||
# corresponding with masked elements are replaced with a sentinel value
|
||||
|
||||
# return without modifying arrays if none have a mask
|
||||
has_mask = False
|
||||
for sample in samples:
|
||||
mask = getattr(sample, 'mask', False)
|
||||
has_mask = has_mask or np.any(mask)
|
||||
if not has_mask:
|
||||
return samples, None # None means there is no sentinel value
|
||||
|
||||
# Choose a sentinel value. We can't use `np.nan`, because sentinel (masked)
|
||||
# values are always omitted, but there are different nan policies.
|
||||
dtype = np.result_type(*samples)
|
||||
dtype = dtype if np.issubdtype(dtype, np.number) else np.float64
|
||||
for i in range(len(samples)):
|
||||
# Things get more complicated if the arrays are of different types.
|
||||
# We could have different sentinel values for each array, but
|
||||
# the purpose of this code is convenience, not efficiency.
|
||||
samples[i] = samples[i].astype(dtype, copy=False)
|
||||
|
||||
inexact = np.issubdtype(dtype, np.inexact)
|
||||
info = np.finfo if inexact else np.iinfo
|
||||
max_possible, min_possible = info(dtype).max, info(dtype).min
|
||||
nextafter = np.nextafter if inexact else (lambda x, _: x - 1)
|
||||
|
||||
sentinel = max_possible
|
||||
# For simplicity, min_possible/np.infs are not candidate sentinel values
|
||||
while sentinel > min_possible:
|
||||
for sample in samples:
|
||||
if np.any(sample == sentinel): # choose a new sentinel value
|
||||
sentinel = nextafter(sentinel, -np.inf)
|
||||
break
|
||||
else: # when sentinel value is OK, break the while loop
|
||||
break
|
||||
else:
|
||||
message = ("This function replaces masked elements with sentinel "
|
||||
"values, but the data contains all distinct values of this "
|
||||
"data type. Consider promoting the dtype to `np.float64`.")
|
||||
raise ValueError(message)
|
||||
|
||||
# replace masked elements with sentinel value
|
||||
out_samples = []
|
||||
for sample in samples:
|
||||
mask = getattr(sample, 'mask', None)
|
||||
if mask is not None: # turn all masked arrays into sentinel arrays
|
||||
mask = np.broadcast_to(mask, sample.shape)
|
||||
sample = sample.data.copy() if np.any(mask) else sample.data
|
||||
sample = np.asarray(sample) # `sample.data` could be a memoryview?
|
||||
sample[mask] = sentinel
|
||||
out_samples.append(sample)
|
||||
|
||||
return out_samples, sentinel
|
||||
|
||||
|
||||
def _check_empty_inputs(samples, axis):
|
||||
"""
|
||||
Check for empty sample; return appropriate output for a vectorized hypotest
|
||||
"""
|
||||
# if none of the samples are empty, we need to perform the test
|
||||
if not any(sample.size == 0 for sample in samples):
|
||||
return None
|
||||
# otherwise, the statistic and p-value will be either empty arrays or
|
||||
# arrays with NaNs. Produce the appropriate array and return it.
|
||||
output_shape = _broadcast_array_shapes_remove_axis(samples, axis)
|
||||
output = np.ones(output_shape) * _get_nan(*samples)
|
||||
return output
|
||||
|
||||
|
||||
def _add_reduced_axes(res, reduced_axes, keepdims):
|
||||
"""
|
||||
Add reduced axes back to all the arrays in the result object
|
||||
if keepdims = True.
|
||||
"""
|
||||
return ([np.expand_dims(output, reduced_axes) for output in res]
|
||||
if keepdims else res)
|
||||
|
||||
|
||||
# Standard docstring / signature entries for `axis`, `nan_policy`, `keepdims`
|
||||
_name = 'axis'
|
||||
_desc = (
|
||||
"""If an int, the axis of the input along which to compute the statistic.
|
||||
The statistic of each axis-slice (e.g. row) of the input will appear in a
|
||||
corresponding element of the output.
|
||||
If ``None``, the input will be raveled before computing the statistic."""
|
||||
.split('\n'))
|
||||
|
||||
|
||||
def _get_axis_params(default_axis=0, _name=_name, _desc=_desc): # bind NOW
|
||||
_type = f"int or None, default: {default_axis}"
|
||||
_axis_parameter_doc = Parameter(_name, _type, _desc)
|
||||
_axis_parameter = inspect.Parameter(_name,
|
||||
inspect.Parameter.KEYWORD_ONLY,
|
||||
default=default_axis)
|
||||
return _axis_parameter_doc, _axis_parameter
|
||||
|
||||
|
||||
_name = 'nan_policy'
|
||||
_type = "{'propagate', 'omit', 'raise'}"
|
||||
_desc = (
|
||||
"""Defines how to handle input NaNs.
|
||||
|
||||
- ``propagate``: if a NaN is present in the axis slice (e.g. row) along
|
||||
which the statistic is computed, the corresponding entry of the output
|
||||
will be NaN.
|
||||
- ``omit``: NaNs will be omitted when performing the calculation.
|
||||
If insufficient data remains in the axis slice along which the
|
||||
statistic is computed, the corresponding entry of the output will be
|
||||
NaN.
|
||||
- ``raise``: if a NaN is present, a ``ValueError`` will be raised."""
|
||||
.split('\n'))
|
||||
_nan_policy_parameter_doc = Parameter(_name, _type, _desc)
|
||||
_nan_policy_parameter = inspect.Parameter(_name,
|
||||
inspect.Parameter.KEYWORD_ONLY,
|
||||
default='propagate')
|
||||
|
||||
_name = 'keepdims'
|
||||
_type = "bool, default: False"
|
||||
_desc = (
|
||||
"""If this is set to True, the axes which are reduced are left
|
||||
in the result as dimensions with size one. With this option,
|
||||
the result will broadcast correctly against the input array."""
|
||||
.split('\n'))
|
||||
_keepdims_parameter_doc = Parameter(_name, _type, _desc)
|
||||
_keepdims_parameter = inspect.Parameter(_name,
|
||||
inspect.Parameter.KEYWORD_ONLY,
|
||||
default=False)
|
||||
|
||||
_standard_note_addition = (
|
||||
"""\nBeginning in SciPy 1.9, ``np.matrix`` inputs (not recommended for new
|
||||
code) are converted to ``np.ndarray`` before the calculation is performed. In
|
||||
this case, the output will be a scalar or ``np.ndarray`` of appropriate shape
|
||||
rather than a 2D ``np.matrix``. Similarly, while masked elements of masked
|
||||
arrays are ignored, the output will be a scalar or ``np.ndarray`` rather than a
|
||||
masked array with ``mask=False``.""").split('\n')
|
||||
|
||||
|
||||
def _axis_nan_policy_factory(tuple_to_result, default_axis=0,
|
||||
n_samples=1, paired=False,
|
||||
result_to_tuple=None, too_small=0,
|
||||
n_outputs=2, kwd_samples=[], override=None):
|
||||
"""Factory for a wrapper that adds axis/nan_policy params to a function.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
tuple_to_result : callable
|
||||
Callable that returns an object of the type returned by the function
|
||||
being wrapped (e.g. the namedtuple or dataclass returned by a
|
||||
statistical test) provided the separate components (e.g. statistic,
|
||||
pvalue).
|
||||
default_axis : int, default: 0
|
||||
The default value of the axis argument. Standard is 0 except when
|
||||
backwards compatibility demands otherwise (e.g. `None`).
|
||||
n_samples : int or callable, default: 1
|
||||
The number of data samples accepted by the function
|
||||
(e.g. `mannwhitneyu`), a callable that accepts a dictionary of
|
||||
parameters passed into the function and returns the number of data
|
||||
samples (e.g. `wilcoxon`), or `None` to indicate an arbitrary number
|
||||
of samples (e.g. `kruskal`).
|
||||
paired : {False, True}
|
||||
Whether the function being wrapped treats the samples as paired (i.e.
|
||||
corresponding elements of each sample should be considered as different
|
||||
components of the same sample.)
|
||||
result_to_tuple : callable, optional
|
||||
Function that unpacks the results of the function being wrapped into
|
||||
a tuple. This is essentially the inverse of `tuple_to_result`. Default
|
||||
is `None`, which is appropriate for statistical tests that return a
|
||||
statistic, pvalue tuple (rather than, e.g., a non-iterable datalass).
|
||||
too_small : int or callable, default: 0
|
||||
The largest unnacceptably small sample for the function being wrapped.
|
||||
For example, some functions require samples of size two or more or they
|
||||
raise an error. This argument prevents the error from being raised when
|
||||
input is not 1D and instead places a NaN in the corresponding element
|
||||
of the result. If callable, it must accept a list of samples, axis,
|
||||
and a dictionary of keyword arguments passed to the wrapper function as
|
||||
arguments and return a bool indicating weather the samples passed are
|
||||
too small.
|
||||
n_outputs : int or callable, default: 2
|
||||
The number of outputs produced by the function given 1d sample(s). For
|
||||
example, hypothesis tests that return a namedtuple or result object
|
||||
with attributes ``statistic`` and ``pvalue`` use the default
|
||||
``n_outputs=2``; summary statistics with scalar output use
|
||||
``n_outputs=1``. Alternatively, may be a callable that accepts a
|
||||
dictionary of arguments passed into the wrapped function and returns
|
||||
the number of outputs corresponding with those arguments.
|
||||
kwd_samples : sequence, default: []
|
||||
The names of keyword parameters that should be treated as samples. For
|
||||
example, `gmean` accepts as its first argument a sample `a` but
|
||||
also `weights` as a fourth, optional keyword argument. In this case, we
|
||||
use `n_samples=1` and kwd_samples=['weights'].
|
||||
override : dict, default: {'vectorization': False, 'nan_propagation': True}
|
||||
Pass a dictionary with ``'vectorization': True`` to ensure that the
|
||||
decorator overrides the function's behavior for multimensional input.
|
||||
Use ``'nan_propagation': False`` to ensure that the decorator does not
|
||||
override the function's behavior for ``nan_policy='propagate'``.
|
||||
(See `scipy.stats.mode`, for example.)
|
||||
"""
|
||||
# Specify which existing behaviors the decorator must override
|
||||
temp = override or {}
|
||||
override = {'vectorization': False,
|
||||
'nan_propagation': True}
|
||||
override.update(temp)
|
||||
|
||||
if result_to_tuple is None:
|
||||
def result_to_tuple(res):
|
||||
return res
|
||||
|
||||
if not callable(too_small):
|
||||
def is_too_small(samples, *ts_args, axis=-1, **ts_kwargs):
|
||||
for sample in samples:
|
||||
if sample.shape[axis] <= too_small:
|
||||
return True
|
||||
return False
|
||||
else:
|
||||
is_too_small = too_small
|
||||
|
||||
def axis_nan_policy_decorator(hypotest_fun_in):
|
||||
@wraps(hypotest_fun_in)
|
||||
def axis_nan_policy_wrapper(*args, _no_deco=False, **kwds):
|
||||
|
||||
if _no_deco: # for testing, decorator does nothing
|
||||
return hypotest_fun_in(*args, **kwds)
|
||||
|
||||
# We need to be flexible about whether position or keyword
|
||||
# arguments are used, but we need to make sure users don't pass
|
||||
# both for the same parameter. To complicate matters, some
|
||||
# functions accept samples with *args, and some functions already
|
||||
# accept `axis` and `nan_policy` as positional arguments.
|
||||
# The strategy is to make sure that there is no duplication
|
||||
# between `args` and `kwds`, combine the two into `kwds`, then
|
||||
# the samples, `nan_policy`, and `axis` from `kwds`, as they are
|
||||
# dealt with separately.
|
||||
|
||||
# Check for intersection between positional and keyword args
|
||||
params = list(inspect.signature(hypotest_fun_in).parameters)
|
||||
if n_samples is None:
|
||||
# Give unique names to each positional sample argument
|
||||
# Note that *args can't be provided as a keyword argument
|
||||
params = [f"arg{i}" for i in range(len(args))] + params[1:]
|
||||
|
||||
# raise if there are too many positional args
|
||||
maxarg = (np.inf if inspect.getfullargspec(hypotest_fun_in).varargs
|
||||
else len(inspect.getfullargspec(hypotest_fun_in).args))
|
||||
if len(args) > maxarg: # let the function raise the right error
|
||||
hypotest_fun_in(*args, **kwds)
|
||||
|
||||
# raise if multiple values passed for same parameter
|
||||
d_args = dict(zip(params, args))
|
||||
intersection = set(d_args) & set(kwds)
|
||||
if intersection: # let the function raise the right error
|
||||
hypotest_fun_in(*args, **kwds)
|
||||
|
||||
# Consolidate other positional and keyword args into `kwds`
|
||||
kwds.update(d_args)
|
||||
|
||||
# rename avoids UnboundLocalError
|
||||
if callable(n_samples):
|
||||
# Future refactoring idea: no need for callable n_samples.
|
||||
# Just replace `n_samples` and `kwd_samples` with a single
|
||||
# list of the names of all samples, and treat all of them
|
||||
# as `kwd_samples` are treated below.
|
||||
n_samp = n_samples(kwds)
|
||||
else:
|
||||
n_samp = n_samples or len(args)
|
||||
|
||||
# get the number of outputs
|
||||
n_out = n_outputs # rename to avoid UnboundLocalError
|
||||
if callable(n_out):
|
||||
n_out = n_out(kwds)
|
||||
|
||||
# If necessary, rearrange function signature: accept other samples
|
||||
# as positional args right after the first n_samp args
|
||||
kwd_samp = [name for name in kwd_samples
|
||||
if kwds.get(name, None) is not None]
|
||||
n_kwd_samp = len(kwd_samp)
|
||||
if not kwd_samp:
|
||||
hypotest_fun_out = hypotest_fun_in
|
||||
else:
|
||||
def hypotest_fun_out(*samples, **kwds):
|
||||
new_kwds = dict(zip(kwd_samp, samples[n_samp:]))
|
||||
kwds.update(new_kwds)
|
||||
return hypotest_fun_in(*samples[:n_samp], **kwds)
|
||||
|
||||
# Extract the things we need here
|
||||
try: # if something is missing
|
||||
samples = [np.atleast_1d(kwds.pop(param))
|
||||
for param in (params[:n_samp] + kwd_samp)]
|
||||
except KeyError: # let the function raise the right error
|
||||
# might need to revisit this if required arg is not a "sample"
|
||||
hypotest_fun_in(*args, **kwds)
|
||||
vectorized = True if 'axis' in params else False
|
||||
vectorized = vectorized and not override['vectorization']
|
||||
axis = kwds.pop('axis', default_axis)
|
||||
nan_policy = kwds.pop('nan_policy', 'propagate')
|
||||
keepdims = kwds.pop("keepdims", False)
|
||||
del args # avoid the possibility of passing both `args` and `kwds`
|
||||
|
||||
# convert masked arrays to regular arrays with sentinel values
|
||||
samples, sentinel = _masked_arrays_2_sentinel_arrays(samples)
|
||||
|
||||
# standardize to always work along last axis
|
||||
reduced_axes = axis
|
||||
if axis is None:
|
||||
if samples:
|
||||
# when axis=None, take the maximum of all dimensions since
|
||||
# all the dimensions are reduced.
|
||||
n_dims = np.max([sample.ndim for sample in samples])
|
||||
reduced_axes = tuple(range(n_dims))
|
||||
samples = [np.asarray(sample.ravel()) for sample in samples]
|
||||
else:
|
||||
samples = _broadcast_arrays(samples, axis=axis)
|
||||
axis = np.atleast_1d(axis)
|
||||
n_axes = len(axis)
|
||||
# move all axes in `axis` to the end to be raveled
|
||||
samples = [np.moveaxis(sample, axis, range(-len(axis), 0))
|
||||
for sample in samples]
|
||||
shapes = [sample.shape for sample in samples]
|
||||
# New shape is unchanged for all axes _not_ in `axis`
|
||||
# At the end, we append the product of the shapes of the axes
|
||||
# in `axis`. Appending -1 doesn't work for zero-size arrays!
|
||||
new_shapes = [shape[:-n_axes] + (np.prod(shape[-n_axes:]),)
|
||||
for shape in shapes]
|
||||
samples = [sample.reshape(new_shape)
|
||||
for sample, new_shape in zip(samples, new_shapes)]
|
||||
axis = -1 # work over the last axis
|
||||
NaN = _get_nan(*samples)
|
||||
|
||||
# if axis is not needed, just handle nan_policy and return
|
||||
ndims = np.array([sample.ndim for sample in samples])
|
||||
if np.all(ndims <= 1):
|
||||
# Addresses nan_policy == "raise"
|
||||
if nan_policy != 'propagate' or override['nan_propagation']:
|
||||
contains_nan = [_contains_nan(sample, nan_policy)[0]
|
||||
for sample in samples]
|
||||
else:
|
||||
# Behave as though there are no NaNs (even if there are)
|
||||
contains_nan = [False]*len(samples)
|
||||
|
||||
# Addresses nan_policy == "propagate"
|
||||
if any(contains_nan) and (nan_policy == 'propagate'
|
||||
and override['nan_propagation']):
|
||||
res = np.full(n_out, NaN)
|
||||
res = _add_reduced_axes(res, reduced_axes, keepdims)
|
||||
return tuple_to_result(*res)
|
||||
|
||||
# Addresses nan_policy == "omit"
|
||||
if any(contains_nan) and nan_policy == 'omit':
|
||||
# consider passing in contains_nan
|
||||
samples = _remove_nans(samples, paired)
|
||||
|
||||
# ideally, this is what the behavior would be:
|
||||
# if is_too_small(samples):
|
||||
# return tuple_to_result(NaN, NaN)
|
||||
# but some existing functions raise exceptions, and changing
|
||||
# behavior of those would break backward compatibility.
|
||||
|
||||
if sentinel:
|
||||
samples = _remove_sentinel(samples, paired, sentinel)
|
||||
res = hypotest_fun_out(*samples, **kwds)
|
||||
res = result_to_tuple(res)
|
||||
res = _add_reduced_axes(res, reduced_axes, keepdims)
|
||||
return tuple_to_result(*res)
|
||||
|
||||
# check for empty input
|
||||
# ideally, move this to the top, but some existing functions raise
|
||||
# exceptions for empty input, so overriding it would break
|
||||
# backward compatibility.
|
||||
empty_output = _check_empty_inputs(samples, axis)
|
||||
# only return empty output if zero sized input is too small.
|
||||
if (
|
||||
empty_output is not None
|
||||
and (is_too_small(samples, kwds) or empty_output.size == 0)
|
||||
):
|
||||
res = [empty_output.copy() for i in range(n_out)]
|
||||
res = _add_reduced_axes(res, reduced_axes, keepdims)
|
||||
return tuple_to_result(*res)
|
||||
|
||||
# otherwise, concatenate all samples along axis, remembering where
|
||||
# each separate sample begins
|
||||
lengths = np.array([sample.shape[axis] for sample in samples])
|
||||
split_indices = np.cumsum(lengths)
|
||||
x = _broadcast_concatenate(samples, axis)
|
||||
|
||||
# Addresses nan_policy == "raise"
|
||||
if nan_policy != 'propagate' or override['nan_propagation']:
|
||||
contains_nan, _ = _contains_nan(x, nan_policy)
|
||||
else:
|
||||
contains_nan = False # behave like there are no NaNs
|
||||
|
||||
if vectorized and not contains_nan and not sentinel:
|
||||
res = hypotest_fun_out(*samples, axis=axis, **kwds)
|
||||
res = result_to_tuple(res)
|
||||
res = _add_reduced_axes(res, reduced_axes, keepdims)
|
||||
return tuple_to_result(*res)
|
||||
|
||||
# Addresses nan_policy == "omit"
|
||||
if contains_nan and nan_policy == 'omit':
|
||||
def hypotest_fun(x):
|
||||
samples = np.split(x, split_indices)[:n_samp+n_kwd_samp]
|
||||
samples = _remove_nans(samples, paired)
|
||||
if sentinel:
|
||||
samples = _remove_sentinel(samples, paired, sentinel)
|
||||
if is_too_small(samples, kwds):
|
||||
return np.full(n_out, NaN)
|
||||
return result_to_tuple(hypotest_fun_out(*samples, **kwds))
|
||||
|
||||
# Addresses nan_policy == "propagate"
|
||||
elif (contains_nan and nan_policy == 'propagate'
|
||||
and override['nan_propagation']):
|
||||
def hypotest_fun(x):
|
||||
if np.isnan(x).any():
|
||||
return np.full(n_out, NaN)
|
||||
|
||||
samples = np.split(x, split_indices)[:n_samp+n_kwd_samp]
|
||||
if sentinel:
|
||||
samples = _remove_sentinel(samples, paired, sentinel)
|
||||
if is_too_small(samples, kwds):
|
||||
return np.full(n_out, NaN)
|
||||
return result_to_tuple(hypotest_fun_out(*samples, **kwds))
|
||||
|
||||
else:
|
||||
def hypotest_fun(x):
|
||||
samples = np.split(x, split_indices)[:n_samp+n_kwd_samp]
|
||||
if sentinel:
|
||||
samples = _remove_sentinel(samples, paired, sentinel)
|
||||
if is_too_small(samples, kwds):
|
||||
return np.full(n_out, NaN)
|
||||
return result_to_tuple(hypotest_fun_out(*samples, **kwds))
|
||||
|
||||
x = np.moveaxis(x, axis, 0)
|
||||
res = np.apply_along_axis(hypotest_fun, axis=0, arr=x)
|
||||
res = _add_reduced_axes(res, reduced_axes, keepdims)
|
||||
return tuple_to_result(*res)
|
||||
|
||||
_axis_parameter_doc, _axis_parameter = _get_axis_params(default_axis)
|
||||
doc = FunctionDoc(axis_nan_policy_wrapper)
|
||||
parameter_names = [param.name for param in doc['Parameters']]
|
||||
if 'axis' in parameter_names:
|
||||
doc['Parameters'][parameter_names.index('axis')] = (
|
||||
_axis_parameter_doc)
|
||||
else:
|
||||
doc['Parameters'].append(_axis_parameter_doc)
|
||||
if 'nan_policy' in parameter_names:
|
||||
doc['Parameters'][parameter_names.index('nan_policy')] = (
|
||||
_nan_policy_parameter_doc)
|
||||
else:
|
||||
doc['Parameters'].append(_nan_policy_parameter_doc)
|
||||
if 'keepdims' in parameter_names:
|
||||
doc['Parameters'][parameter_names.index('keepdims')] = (
|
||||
_keepdims_parameter_doc)
|
||||
else:
|
||||
doc['Parameters'].append(_keepdims_parameter_doc)
|
||||
doc['Notes'] += _standard_note_addition
|
||||
doc = str(doc).split("\n", 1)[1] # remove signature
|
||||
axis_nan_policy_wrapper.__doc__ = str(doc)
|
||||
|
||||
sig = inspect.signature(axis_nan_policy_wrapper)
|
||||
parameters = sig.parameters
|
||||
parameter_list = list(parameters.values())
|
||||
if 'axis' not in parameters:
|
||||
parameter_list.append(_axis_parameter)
|
||||
if 'nan_policy' not in parameters:
|
||||
parameter_list.append(_nan_policy_parameter)
|
||||
if 'keepdims' not in parameters:
|
||||
parameter_list.append(_keepdims_parameter)
|
||||
sig = sig.replace(parameters=parameter_list)
|
||||
axis_nan_policy_wrapper.__signature__ = sig
|
||||
|
||||
return axis_nan_policy_wrapper
|
||||
return axis_nan_policy_decorator
|
||||
Binary file not shown.
Binary file not shown.
@@ -0,0 +1,27 @@
|
||||
# Declare the class with cdef
|
||||
cdef extern from "biasedurn/stocc.h" nogil:
|
||||
cdef cppclass CFishersNCHypergeometric:
|
||||
CFishersNCHypergeometric(int, int, int, double, double) except +
|
||||
int mode()
|
||||
double mean()
|
||||
double variance()
|
||||
double probability(int x)
|
||||
double moments(double * mean, double * var)
|
||||
|
||||
cdef cppclass CWalleniusNCHypergeometric:
|
||||
CWalleniusNCHypergeometric() except +
|
||||
CWalleniusNCHypergeometric(int, int, int, double, double) except +
|
||||
int mode()
|
||||
double mean()
|
||||
double variance()
|
||||
double probability(int x)
|
||||
double moments(double * mean, double * var)
|
||||
|
||||
cdef cppclass StochasticLib3:
|
||||
StochasticLib3(int seed) except +
|
||||
double Random() except +
|
||||
void SetAccuracy(double accur)
|
||||
int FishersNCHyp (int n, int m, int N, double odds) except +
|
||||
int WalleniusNCHyp (int n, int m, int N, double odds) except +
|
||||
double(*next_double)()
|
||||
double(*next_normal)(const double m, const double s)
|
||||
@@ -0,0 +1,795 @@
|
||||
import builtins
|
||||
from warnings import catch_warnings, simplefilter
|
||||
import numpy as np
|
||||
from operator import index
|
||||
from collections import namedtuple
|
||||
|
||||
__all__ = ['binned_statistic',
|
||||
'binned_statistic_2d',
|
||||
'binned_statistic_dd']
|
||||
|
||||
|
||||
BinnedStatisticResult = namedtuple('BinnedStatisticResult',
|
||||
('statistic', 'bin_edges', 'binnumber'))
|
||||
|
||||
|
||||
def binned_statistic(x, values, statistic='mean',
|
||||
bins=10, range=None):
|
||||
"""
|
||||
Compute a binned statistic for one or more sets of data.
|
||||
|
||||
This is a generalization of a histogram function. A histogram divides
|
||||
the space into bins, and returns the count of the number of points in
|
||||
each bin. This function allows the computation of the sum, mean, median,
|
||||
or other statistic of the values (or set of values) within each bin.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x : (N,) array_like
|
||||
A sequence of values to be binned.
|
||||
values : (N,) array_like or list of (N,) array_like
|
||||
The data on which the statistic will be computed. This must be
|
||||
the same shape as `x`, or a set of sequences - each the same shape as
|
||||
`x`. If `values` is a set of sequences, the statistic will be computed
|
||||
on each independently.
|
||||
statistic : string or callable, optional
|
||||
The statistic to compute (default is 'mean').
|
||||
The following statistics are available:
|
||||
|
||||
* 'mean' : compute the mean of values for points within each bin.
|
||||
Empty bins will be represented by NaN.
|
||||
* 'std' : compute the standard deviation within each bin. This
|
||||
is implicitly calculated with ddof=0.
|
||||
* 'median' : compute the median of values for points within each
|
||||
bin. Empty bins will be represented by NaN.
|
||||
* 'count' : compute the count of points within each bin. This is
|
||||
identical to an unweighted histogram. `values` array is not
|
||||
referenced.
|
||||
* 'sum' : compute the sum of values for points within each bin.
|
||||
This is identical to a weighted histogram.
|
||||
* 'min' : compute the minimum of values for points within each bin.
|
||||
Empty bins will be represented by NaN.
|
||||
* 'max' : compute the maximum of values for point within each bin.
|
||||
Empty bins will be represented by NaN.
|
||||
* function : a user-defined function which takes a 1D array of
|
||||
values, and outputs a single numerical statistic. This function
|
||||
will be called on the values in each bin. Empty bins will be
|
||||
represented by function([]), or NaN if this returns an error.
|
||||
|
||||
bins : int or sequence of scalars, optional
|
||||
If `bins` is an int, it defines the number of equal-width bins in the
|
||||
given range (10 by default). If `bins` is a sequence, it defines the
|
||||
bin edges, including the rightmost edge, allowing for non-uniform bin
|
||||
widths. Values in `x` that are smaller than lowest bin edge are
|
||||
assigned to bin number 0, values beyond the highest bin are assigned to
|
||||
``bins[-1]``. If the bin edges are specified, the number of bins will
|
||||
be, (nx = len(bins)-1).
|
||||
range : (float, float) or [(float, float)], optional
|
||||
The lower and upper range of the bins. If not provided, range
|
||||
is simply ``(x.min(), x.max())``. Values outside the range are
|
||||
ignored.
|
||||
|
||||
Returns
|
||||
-------
|
||||
statistic : array
|
||||
The values of the selected statistic in each bin.
|
||||
bin_edges : array of dtype float
|
||||
Return the bin edges ``(length(statistic)+1)``.
|
||||
binnumber: 1-D ndarray of ints
|
||||
Indices of the bins (corresponding to `bin_edges`) in which each value
|
||||
of `x` belongs. Same length as `values`. A binnumber of `i` means the
|
||||
corresponding value is between (bin_edges[i-1], bin_edges[i]).
|
||||
|
||||
See Also
|
||||
--------
|
||||
numpy.digitize, numpy.histogram, binned_statistic_2d, binned_statistic_dd
|
||||
|
||||
Notes
|
||||
-----
|
||||
All but the last (righthand-most) bin is half-open. In other words, if
|
||||
`bins` is ``[1, 2, 3, 4]``, then the first bin is ``[1, 2)`` (including 1,
|
||||
but excluding 2) and the second ``[2, 3)``. The last bin, however, is
|
||||
``[3, 4]``, which *includes* 4.
|
||||
|
||||
.. versionadded:: 0.11.0
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> import numpy as np
|
||||
>>> from scipy import stats
|
||||
>>> import matplotlib.pyplot as plt
|
||||
|
||||
First some basic examples:
|
||||
|
||||
Create two evenly spaced bins in the range of the given sample, and sum the
|
||||
corresponding values in each of those bins:
|
||||
|
||||
>>> values = [1.0, 1.0, 2.0, 1.5, 3.0]
|
||||
>>> stats.binned_statistic([1, 1, 2, 5, 7], values, 'sum', bins=2)
|
||||
BinnedStatisticResult(statistic=array([4. , 4.5]),
|
||||
bin_edges=array([1., 4., 7.]), binnumber=array([1, 1, 1, 2, 2]))
|
||||
|
||||
Multiple arrays of values can also be passed. The statistic is calculated
|
||||
on each set independently:
|
||||
|
||||
>>> values = [[1.0, 1.0, 2.0, 1.5, 3.0], [2.0, 2.0, 4.0, 3.0, 6.0]]
|
||||
>>> stats.binned_statistic([1, 1, 2, 5, 7], values, 'sum', bins=2)
|
||||
BinnedStatisticResult(statistic=array([[4. , 4.5],
|
||||
[8. , 9. ]]), bin_edges=array([1., 4., 7.]),
|
||||
binnumber=array([1, 1, 1, 2, 2]))
|
||||
|
||||
>>> stats.binned_statistic([1, 2, 1, 2, 4], np.arange(5), statistic='mean',
|
||||
... bins=3)
|
||||
BinnedStatisticResult(statistic=array([1., 2., 4.]),
|
||||
bin_edges=array([1., 2., 3., 4.]),
|
||||
binnumber=array([1, 2, 1, 2, 3]))
|
||||
|
||||
As a second example, we now generate some random data of sailing boat speed
|
||||
as a function of wind speed, and then determine how fast our boat is for
|
||||
certain wind speeds:
|
||||
|
||||
>>> rng = np.random.default_rng()
|
||||
>>> windspeed = 8 * rng.random(500)
|
||||
>>> boatspeed = .3 * windspeed**.5 + .2 * rng.random(500)
|
||||
>>> bin_means, bin_edges, binnumber = stats.binned_statistic(windspeed,
|
||||
... boatspeed, statistic='median', bins=[1,2,3,4,5,6,7])
|
||||
>>> plt.figure()
|
||||
>>> plt.plot(windspeed, boatspeed, 'b.', label='raw data')
|
||||
>>> plt.hlines(bin_means, bin_edges[:-1], bin_edges[1:], colors='g', lw=5,
|
||||
... label='binned statistic of data')
|
||||
>>> plt.legend()
|
||||
|
||||
Now we can use ``binnumber`` to select all datapoints with a windspeed
|
||||
below 1:
|
||||
|
||||
>>> low_boatspeed = boatspeed[binnumber == 0]
|
||||
|
||||
As a final example, we will use ``bin_edges`` and ``binnumber`` to make a
|
||||
plot of a distribution that shows the mean and distribution around that
|
||||
mean per bin, on top of a regular histogram and the probability
|
||||
distribution function:
|
||||
|
||||
>>> x = np.linspace(0, 5, num=500)
|
||||
>>> x_pdf = stats.maxwell.pdf(x)
|
||||
>>> samples = stats.maxwell.rvs(size=10000)
|
||||
|
||||
>>> bin_means, bin_edges, binnumber = stats.binned_statistic(x, x_pdf,
|
||||
... statistic='mean', bins=25)
|
||||
>>> bin_width = (bin_edges[1] - bin_edges[0])
|
||||
>>> bin_centers = bin_edges[1:] - bin_width/2
|
||||
|
||||
>>> plt.figure()
|
||||
>>> plt.hist(samples, bins=50, density=True, histtype='stepfilled',
|
||||
... alpha=0.2, label='histogram of data')
|
||||
>>> plt.plot(x, x_pdf, 'r-', label='analytical pdf')
|
||||
>>> plt.hlines(bin_means, bin_edges[:-1], bin_edges[1:], colors='g', lw=2,
|
||||
... label='binned statistic of data')
|
||||
>>> plt.plot((binnumber - 0.5) * bin_width, x_pdf, 'g.', alpha=0.5)
|
||||
>>> plt.legend(fontsize=10)
|
||||
>>> plt.show()
|
||||
|
||||
"""
|
||||
try:
|
||||
N = len(bins)
|
||||
except TypeError:
|
||||
N = 1
|
||||
|
||||
if N != 1:
|
||||
bins = [np.asarray(bins, float)]
|
||||
|
||||
if range is not None:
|
||||
if len(range) == 2:
|
||||
range = [range]
|
||||
|
||||
medians, edges, binnumbers = binned_statistic_dd(
|
||||
[x], values, statistic, bins, range)
|
||||
|
||||
return BinnedStatisticResult(medians, edges[0], binnumbers)
|
||||
|
||||
|
||||
BinnedStatistic2dResult = namedtuple('BinnedStatistic2dResult',
|
||||
('statistic', 'x_edge', 'y_edge',
|
||||
'binnumber'))
|
||||
|
||||
|
||||
def binned_statistic_2d(x, y, values, statistic='mean',
|
||||
bins=10, range=None, expand_binnumbers=False):
|
||||
"""
|
||||
Compute a bidimensional binned statistic for one or more sets of data.
|
||||
|
||||
This is a generalization of a histogram2d function. A histogram divides
|
||||
the space into bins, and returns the count of the number of points in
|
||||
each bin. This function allows the computation of the sum, mean, median,
|
||||
or other statistic of the values (or set of values) within each bin.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x : (N,) array_like
|
||||
A sequence of values to be binned along the first dimension.
|
||||
y : (N,) array_like
|
||||
A sequence of values to be binned along the second dimension.
|
||||
values : (N,) array_like or list of (N,) array_like
|
||||
The data on which the statistic will be computed. This must be
|
||||
the same shape as `x`, or a list of sequences - each with the same
|
||||
shape as `x`. If `values` is such a list, the statistic will be
|
||||
computed on each independently.
|
||||
statistic : string or callable, optional
|
||||
The statistic to compute (default is 'mean').
|
||||
The following statistics are available:
|
||||
|
||||
* 'mean' : compute the mean of values for points within each bin.
|
||||
Empty bins will be represented by NaN.
|
||||
* 'std' : compute the standard deviation within each bin. This
|
||||
is implicitly calculated with ddof=0.
|
||||
* 'median' : compute the median of values for points within each
|
||||
bin. Empty bins will be represented by NaN.
|
||||
* 'count' : compute the count of points within each bin. This is
|
||||
identical to an unweighted histogram. `values` array is not
|
||||
referenced.
|
||||
* 'sum' : compute the sum of values for points within each bin.
|
||||
This is identical to a weighted histogram.
|
||||
* 'min' : compute the minimum of values for points within each bin.
|
||||
Empty bins will be represented by NaN.
|
||||
* 'max' : compute the maximum of values for point within each bin.
|
||||
Empty bins will be represented by NaN.
|
||||
* function : a user-defined function which takes a 1D array of
|
||||
values, and outputs a single numerical statistic. This function
|
||||
will be called on the values in each bin. Empty bins will be
|
||||
represented by function([]), or NaN if this returns an error.
|
||||
|
||||
bins : int or [int, int] or array_like or [array, array], optional
|
||||
The bin specification:
|
||||
|
||||
* the number of bins for the two dimensions (nx = ny = bins),
|
||||
* the number of bins in each dimension (nx, ny = bins),
|
||||
* the bin edges for the two dimensions (x_edge = y_edge = bins),
|
||||
* the bin edges in each dimension (x_edge, y_edge = bins).
|
||||
|
||||
If the bin edges are specified, the number of bins will be,
|
||||
(nx = len(x_edge)-1, ny = len(y_edge)-1).
|
||||
|
||||
range : (2,2) array_like, optional
|
||||
The leftmost and rightmost edges of the bins along each dimension
|
||||
(if not specified explicitly in the `bins` parameters):
|
||||
[[xmin, xmax], [ymin, ymax]]. All values outside of this range will be
|
||||
considered outliers and not tallied in the histogram.
|
||||
expand_binnumbers : bool, optional
|
||||
'False' (default): the returned `binnumber` is a shape (N,) array of
|
||||
linearized bin indices.
|
||||
'True': the returned `binnumber` is 'unraveled' into a shape (2,N)
|
||||
ndarray, where each row gives the bin numbers in the corresponding
|
||||
dimension.
|
||||
See the `binnumber` returned value, and the `Examples` section.
|
||||
|
||||
.. versionadded:: 0.17.0
|
||||
|
||||
Returns
|
||||
-------
|
||||
statistic : (nx, ny) ndarray
|
||||
The values of the selected statistic in each two-dimensional bin.
|
||||
x_edge : (nx + 1) ndarray
|
||||
The bin edges along the first dimension.
|
||||
y_edge : (ny + 1) ndarray
|
||||
The bin edges along the second dimension.
|
||||
binnumber : (N,) array of ints or (2,N) ndarray of ints
|
||||
This assigns to each element of `sample` an integer that represents the
|
||||
bin in which this observation falls. The representation depends on the
|
||||
`expand_binnumbers` argument. See `Notes` for details.
|
||||
|
||||
|
||||
See Also
|
||||
--------
|
||||
numpy.digitize, numpy.histogram2d, binned_statistic, binned_statistic_dd
|
||||
|
||||
Notes
|
||||
-----
|
||||
Binedges:
|
||||
All but the last (righthand-most) bin is half-open. In other words, if
|
||||
`bins` is ``[1, 2, 3, 4]``, then the first bin is ``[1, 2)`` (including 1,
|
||||
but excluding 2) and the second ``[2, 3)``. The last bin, however, is
|
||||
``[3, 4]``, which *includes* 4.
|
||||
|
||||
`binnumber`:
|
||||
This returned argument assigns to each element of `sample` an integer that
|
||||
represents the bin in which it belongs. The representation depends on the
|
||||
`expand_binnumbers` argument. If 'False' (default): The returned
|
||||
`binnumber` is a shape (N,) array of linearized indices mapping each
|
||||
element of `sample` to its corresponding bin (using row-major ordering).
|
||||
Note that the returned linearized bin indices are used for an array with
|
||||
extra bins on the outer binedges to capture values outside of the defined
|
||||
bin bounds.
|
||||
If 'True': The returned `binnumber` is a shape (2,N) ndarray where
|
||||
each row indicates bin placements for each dimension respectively. In each
|
||||
dimension, a binnumber of `i` means the corresponding value is between
|
||||
(D_edge[i-1], D_edge[i]), where 'D' is either 'x' or 'y'.
|
||||
|
||||
.. versionadded:: 0.11.0
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> from scipy import stats
|
||||
|
||||
Calculate the counts with explicit bin-edges:
|
||||
|
||||
>>> x = [0.1, 0.1, 0.1, 0.6]
|
||||
>>> y = [2.1, 2.6, 2.1, 2.1]
|
||||
>>> binx = [0.0, 0.5, 1.0]
|
||||
>>> biny = [2.0, 2.5, 3.0]
|
||||
>>> ret = stats.binned_statistic_2d(x, y, None, 'count', bins=[binx, biny])
|
||||
>>> ret.statistic
|
||||
array([[2., 1.],
|
||||
[1., 0.]])
|
||||
|
||||
The bin in which each sample is placed is given by the `binnumber`
|
||||
returned parameter. By default, these are the linearized bin indices:
|
||||
|
||||
>>> ret.binnumber
|
||||
array([5, 6, 5, 9])
|
||||
|
||||
The bin indices can also be expanded into separate entries for each
|
||||
dimension using the `expand_binnumbers` parameter:
|
||||
|
||||
>>> ret = stats.binned_statistic_2d(x, y, None, 'count', bins=[binx, biny],
|
||||
... expand_binnumbers=True)
|
||||
>>> ret.binnumber
|
||||
array([[1, 1, 1, 2],
|
||||
[1, 2, 1, 1]])
|
||||
|
||||
Which shows that the first three elements belong in the xbin 1, and the
|
||||
fourth into xbin 2; and so on for y.
|
||||
|
||||
"""
|
||||
|
||||
# This code is based on np.histogram2d
|
||||
try:
|
||||
N = len(bins)
|
||||
except TypeError:
|
||||
N = 1
|
||||
|
||||
if N != 1 and N != 2:
|
||||
xedges = yedges = np.asarray(bins, float)
|
||||
bins = [xedges, yedges]
|
||||
|
||||
medians, edges, binnumbers = binned_statistic_dd(
|
||||
[x, y], values, statistic, bins, range,
|
||||
expand_binnumbers=expand_binnumbers)
|
||||
|
||||
return BinnedStatistic2dResult(medians, edges[0], edges[1], binnumbers)
|
||||
|
||||
|
||||
BinnedStatisticddResult = namedtuple('BinnedStatisticddResult',
|
||||
('statistic', 'bin_edges',
|
||||
'binnumber'))
|
||||
|
||||
|
||||
def _bincount(x, weights):
|
||||
if np.iscomplexobj(weights):
|
||||
a = np.bincount(x, np.real(weights))
|
||||
b = np.bincount(x, np.imag(weights))
|
||||
z = a + b*1j
|
||||
|
||||
else:
|
||||
z = np.bincount(x, weights)
|
||||
return z
|
||||
|
||||
|
||||
def binned_statistic_dd(sample, values, statistic='mean',
|
||||
bins=10, range=None, expand_binnumbers=False,
|
||||
binned_statistic_result=None):
|
||||
"""
|
||||
Compute a multidimensional binned statistic for a set of data.
|
||||
|
||||
This is a generalization of a histogramdd function. A histogram divides
|
||||
the space into bins, and returns the count of the number of points in
|
||||
each bin. This function allows the computation of the sum, mean, median,
|
||||
or other statistic of the values within each bin.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
sample : array_like
|
||||
Data to histogram passed as a sequence of N arrays of length D, or
|
||||
as an (N,D) array.
|
||||
values : (N,) array_like or list of (N,) array_like
|
||||
The data on which the statistic will be computed. This must be
|
||||
the same shape as `sample`, or a list of sequences - each with the
|
||||
same shape as `sample`. If `values` is such a list, the statistic
|
||||
will be computed on each independently.
|
||||
statistic : string or callable, optional
|
||||
The statistic to compute (default is 'mean').
|
||||
The following statistics are available:
|
||||
|
||||
* 'mean' : compute the mean of values for points within each bin.
|
||||
Empty bins will be represented by NaN.
|
||||
* 'median' : compute the median of values for points within each
|
||||
bin. Empty bins will be represented by NaN.
|
||||
* 'count' : compute the count of points within each bin. This is
|
||||
identical to an unweighted histogram. `values` array is not
|
||||
referenced.
|
||||
* 'sum' : compute the sum of values for points within each bin.
|
||||
This is identical to a weighted histogram.
|
||||
* 'std' : compute the standard deviation within each bin. This
|
||||
is implicitly calculated with ddof=0. If the number of values
|
||||
within a given bin is 0 or 1, the computed standard deviation value
|
||||
will be 0 for the bin.
|
||||
* 'min' : compute the minimum of values for points within each bin.
|
||||
Empty bins will be represented by NaN.
|
||||
* 'max' : compute the maximum of values for point within each bin.
|
||||
Empty bins will be represented by NaN.
|
||||
* function : a user-defined function which takes a 1D array of
|
||||
values, and outputs a single numerical statistic. This function
|
||||
will be called on the values in each bin. Empty bins will be
|
||||
represented by function([]), or NaN if this returns an error.
|
||||
|
||||
bins : sequence or positive int, optional
|
||||
The bin specification must be in one of the following forms:
|
||||
|
||||
* A sequence of arrays describing the bin edges along each dimension.
|
||||
* The number of bins for each dimension (nx, ny, ... = bins).
|
||||
* The number of bins for all dimensions (nx = ny = ... = bins).
|
||||
range : sequence, optional
|
||||
A sequence of lower and upper bin edges to be used if the edges are
|
||||
not given explicitly in `bins`. Defaults to the minimum and maximum
|
||||
values along each dimension.
|
||||
expand_binnumbers : bool, optional
|
||||
'False' (default): the returned `binnumber` is a shape (N,) array of
|
||||
linearized bin indices.
|
||||
'True': the returned `binnumber` is 'unraveled' into a shape (D,N)
|
||||
ndarray, where each row gives the bin numbers in the corresponding
|
||||
dimension.
|
||||
See the `binnumber` returned value, and the `Examples` section of
|
||||
`binned_statistic_2d`.
|
||||
binned_statistic_result : binnedStatisticddResult
|
||||
Result of a previous call to the function in order to reuse bin edges
|
||||
and bin numbers with new values and/or a different statistic.
|
||||
To reuse bin numbers, `expand_binnumbers` must have been set to False
|
||||
(the default)
|
||||
|
||||
.. versionadded:: 0.17.0
|
||||
|
||||
Returns
|
||||
-------
|
||||
statistic : ndarray, shape(nx1, nx2, nx3,...)
|
||||
The values of the selected statistic in each two-dimensional bin.
|
||||
bin_edges : list of ndarrays
|
||||
A list of D arrays describing the (nxi + 1) bin edges for each
|
||||
dimension.
|
||||
binnumber : (N,) array of ints or (D,N) ndarray of ints
|
||||
This assigns to each element of `sample` an integer that represents the
|
||||
bin in which this observation falls. The representation depends on the
|
||||
`expand_binnumbers` argument. See `Notes` for details.
|
||||
|
||||
|
||||
See Also
|
||||
--------
|
||||
numpy.digitize, numpy.histogramdd, binned_statistic, binned_statistic_2d
|
||||
|
||||
Notes
|
||||
-----
|
||||
Binedges:
|
||||
All but the last (righthand-most) bin is half-open in each dimension. In
|
||||
other words, if `bins` is ``[1, 2, 3, 4]``, then the first bin is
|
||||
``[1, 2)`` (including 1, but excluding 2) and the second ``[2, 3)``. The
|
||||
last bin, however, is ``[3, 4]``, which *includes* 4.
|
||||
|
||||
`binnumber`:
|
||||
This returned argument assigns to each element of `sample` an integer that
|
||||
represents the bin in which it belongs. The representation depends on the
|
||||
`expand_binnumbers` argument. If 'False' (default): The returned
|
||||
`binnumber` is a shape (N,) array of linearized indices mapping each
|
||||
element of `sample` to its corresponding bin (using row-major ordering).
|
||||
If 'True': The returned `binnumber` is a shape (D,N) ndarray where
|
||||
each row indicates bin placements for each dimension respectively. In each
|
||||
dimension, a binnumber of `i` means the corresponding value is between
|
||||
(bin_edges[D][i-1], bin_edges[D][i]), for each dimension 'D'.
|
||||
|
||||
.. versionadded:: 0.11.0
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> import numpy as np
|
||||
>>> from scipy import stats
|
||||
>>> import matplotlib.pyplot as plt
|
||||
>>> from mpl_toolkits.mplot3d import Axes3D
|
||||
|
||||
Take an array of 600 (x, y) coordinates as an example.
|
||||
`binned_statistic_dd` can handle arrays of higher dimension `D`. But a plot
|
||||
of dimension `D+1` is required.
|
||||
|
||||
>>> mu = np.array([0., 1.])
|
||||
>>> sigma = np.array([[1., -0.5],[-0.5, 1.5]])
|
||||
>>> multinormal = stats.multivariate_normal(mu, sigma)
|
||||
>>> data = multinormal.rvs(size=600, random_state=235412)
|
||||
>>> data.shape
|
||||
(600, 2)
|
||||
|
||||
Create bins and count how many arrays fall in each bin:
|
||||
|
||||
>>> N = 60
|
||||
>>> x = np.linspace(-3, 3, N)
|
||||
>>> y = np.linspace(-3, 4, N)
|
||||
>>> ret = stats.binned_statistic_dd(data, np.arange(600), bins=[x, y],
|
||||
... statistic='count')
|
||||
>>> bincounts = ret.statistic
|
||||
|
||||
Set the volume and the location of bars:
|
||||
|
||||
>>> dx = x[1] - x[0]
|
||||
>>> dy = y[1] - y[0]
|
||||
>>> x, y = np.meshgrid(x[:-1]+dx/2, y[:-1]+dy/2)
|
||||
>>> z = 0
|
||||
|
||||
>>> bincounts = bincounts.ravel()
|
||||
>>> x = x.ravel()
|
||||
>>> y = y.ravel()
|
||||
|
||||
>>> fig = plt.figure()
|
||||
>>> ax = fig.add_subplot(111, projection='3d')
|
||||
>>> with np.errstate(divide='ignore'): # silence random axes3d warning
|
||||
... ax.bar3d(x, y, z, dx, dy, bincounts)
|
||||
|
||||
Reuse bin numbers and bin edges with new values:
|
||||
|
||||
>>> ret2 = stats.binned_statistic_dd(data, -np.arange(600),
|
||||
... binned_statistic_result=ret,
|
||||
... statistic='mean')
|
||||
"""
|
||||
known_stats = ['mean', 'median', 'count', 'sum', 'std', 'min', 'max']
|
||||
if not callable(statistic) and statistic not in known_stats:
|
||||
raise ValueError(f'invalid statistic {statistic!r}')
|
||||
|
||||
try:
|
||||
bins = index(bins)
|
||||
except TypeError:
|
||||
# bins is not an integer
|
||||
pass
|
||||
# If bins was an integer-like object, now it is an actual Python int.
|
||||
|
||||
# NOTE: for _bin_edges(), see e.g. gh-11365
|
||||
if isinstance(bins, int) and not np.isfinite(sample).all():
|
||||
raise ValueError(f'{sample!r} contains non-finite values.')
|
||||
|
||||
# `Ndim` is the number of dimensions (e.g. `2` for `binned_statistic_2d`)
|
||||
# `Dlen` is the length of elements along each dimension.
|
||||
# This code is based on np.histogramdd
|
||||
try:
|
||||
# `sample` is an ND-array.
|
||||
Dlen, Ndim = sample.shape
|
||||
except (AttributeError, ValueError):
|
||||
# `sample` is a sequence of 1D arrays.
|
||||
sample = np.atleast_2d(sample).T
|
||||
Dlen, Ndim = sample.shape
|
||||
|
||||
# Store initial shape of `values` to preserve it in the output
|
||||
values = np.asarray(values)
|
||||
input_shape = list(values.shape)
|
||||
# Make sure that `values` is 2D to iterate over rows
|
||||
values = np.atleast_2d(values)
|
||||
Vdim, Vlen = values.shape
|
||||
|
||||
# Make sure `values` match `sample`
|
||||
if statistic != 'count' and Vlen != Dlen:
|
||||
raise AttributeError('The number of `values` elements must match the '
|
||||
'length of each `sample` dimension.')
|
||||
|
||||
try:
|
||||
M = len(bins)
|
||||
if M != Ndim:
|
||||
raise AttributeError('The dimension of bins must be equal '
|
||||
'to the dimension of the sample x.')
|
||||
except TypeError:
|
||||
bins = Ndim * [bins]
|
||||
|
||||
if binned_statistic_result is None:
|
||||
nbin, edges, dedges = _bin_edges(sample, bins, range)
|
||||
binnumbers = _bin_numbers(sample, nbin, edges, dedges)
|
||||
else:
|
||||
edges = binned_statistic_result.bin_edges
|
||||
nbin = np.array([len(edges[i]) + 1 for i in builtins.range(Ndim)])
|
||||
# +1 for outlier bins
|
||||
dedges = [np.diff(edges[i]) for i in builtins.range(Ndim)]
|
||||
binnumbers = binned_statistic_result.binnumber
|
||||
|
||||
# Avoid overflow with double precision. Complex `values` -> `complex128`.
|
||||
result_type = np.result_type(values, np.float64)
|
||||
result = np.empty([Vdim, nbin.prod()], dtype=result_type)
|
||||
|
||||
if statistic in {'mean', np.mean}:
|
||||
result.fill(np.nan)
|
||||
flatcount = _bincount(binnumbers, None)
|
||||
a = flatcount.nonzero()
|
||||
for vv in builtins.range(Vdim):
|
||||
flatsum = _bincount(binnumbers, values[vv])
|
||||
result[vv, a] = flatsum[a] / flatcount[a]
|
||||
elif statistic in {'std', np.std}:
|
||||
result.fill(np.nan)
|
||||
flatcount = _bincount(binnumbers, None)
|
||||
a = flatcount.nonzero()
|
||||
for vv in builtins.range(Vdim):
|
||||
flatsum = _bincount(binnumbers, values[vv])
|
||||
delta = values[vv] - flatsum[binnumbers] / flatcount[binnumbers]
|
||||
std = np.sqrt(
|
||||
_bincount(binnumbers, delta*np.conj(delta))[a] / flatcount[a]
|
||||
)
|
||||
result[vv, a] = std
|
||||
result = np.real(result)
|
||||
elif statistic == 'count':
|
||||
result = np.empty([Vdim, nbin.prod()], dtype=np.float64)
|
||||
result.fill(0)
|
||||
flatcount = _bincount(binnumbers, None)
|
||||
a = np.arange(len(flatcount))
|
||||
result[:, a] = flatcount[np.newaxis, :]
|
||||
elif statistic in {'sum', np.sum}:
|
||||
result.fill(0)
|
||||
for vv in builtins.range(Vdim):
|
||||
flatsum = _bincount(binnumbers, values[vv])
|
||||
a = np.arange(len(flatsum))
|
||||
result[vv, a] = flatsum
|
||||
elif statistic in {'median', np.median}:
|
||||
result.fill(np.nan)
|
||||
for vv in builtins.range(Vdim):
|
||||
i = np.lexsort((values[vv], binnumbers))
|
||||
_, j, counts = np.unique(binnumbers[i],
|
||||
return_index=True, return_counts=True)
|
||||
mid = j + (counts - 1) / 2
|
||||
mid_a = values[vv, i][np.floor(mid).astype(int)]
|
||||
mid_b = values[vv, i][np.ceil(mid).astype(int)]
|
||||
medians = (mid_a + mid_b) / 2
|
||||
result[vv, binnumbers[i][j]] = medians
|
||||
elif statistic in {'min', np.min}:
|
||||
result.fill(np.nan)
|
||||
for vv in builtins.range(Vdim):
|
||||
i = np.argsort(values[vv])[::-1] # Reversed so the min is last
|
||||
result[vv, binnumbers[i]] = values[vv, i]
|
||||
elif statistic in {'max', np.max}:
|
||||
result.fill(np.nan)
|
||||
for vv in builtins.range(Vdim):
|
||||
i = np.argsort(values[vv])
|
||||
result[vv, binnumbers[i]] = values[vv, i]
|
||||
elif callable(statistic):
|
||||
with np.errstate(invalid='ignore'), catch_warnings():
|
||||
simplefilter("ignore", RuntimeWarning)
|
||||
try:
|
||||
null = statistic([])
|
||||
except Exception:
|
||||
null = np.nan
|
||||
if np.iscomplexobj(null):
|
||||
result = result.astype(np.complex128)
|
||||
result.fill(null)
|
||||
try:
|
||||
_calc_binned_statistic(
|
||||
Vdim, binnumbers, result, values, statistic
|
||||
)
|
||||
except ValueError:
|
||||
result = result.astype(np.complex128)
|
||||
_calc_binned_statistic(
|
||||
Vdim, binnumbers, result, values, statistic
|
||||
)
|
||||
|
||||
# Shape into a proper matrix
|
||||
result = result.reshape(np.append(Vdim, nbin))
|
||||
|
||||
# Remove outliers (indices 0 and -1 for each bin-dimension).
|
||||
core = tuple([slice(None)] + Ndim * [slice(1, -1)])
|
||||
result = result[core]
|
||||
|
||||
# Unravel binnumbers into an ndarray, each row the bins for each dimension
|
||||
if expand_binnumbers and Ndim > 1:
|
||||
binnumbers = np.asarray(np.unravel_index(binnumbers, nbin))
|
||||
|
||||
if np.any(result.shape[1:] != nbin - 2):
|
||||
raise RuntimeError('Internal Shape Error')
|
||||
|
||||
# Reshape to have output (`result`) match input (`values`) shape
|
||||
result = result.reshape(input_shape[:-1] + list(nbin-2))
|
||||
|
||||
return BinnedStatisticddResult(result, edges, binnumbers)
|
||||
|
||||
|
||||
def _calc_binned_statistic(Vdim, bin_numbers, result, values, stat_func):
|
||||
unique_bin_numbers = np.unique(bin_numbers)
|
||||
for vv in builtins.range(Vdim):
|
||||
bin_map = _create_binned_data(bin_numbers, unique_bin_numbers,
|
||||
values, vv)
|
||||
for i in unique_bin_numbers:
|
||||
stat = stat_func(np.array(bin_map[i]))
|
||||
if np.iscomplexobj(stat) and not np.iscomplexobj(result):
|
||||
raise ValueError("The statistic function returns complex ")
|
||||
result[vv, i] = stat
|
||||
|
||||
|
||||
def _create_binned_data(bin_numbers, unique_bin_numbers, values, vv):
|
||||
""" Create hashmap of bin ids to values in bins
|
||||
key: bin number
|
||||
value: list of binned data
|
||||
"""
|
||||
bin_map = dict()
|
||||
for i in unique_bin_numbers:
|
||||
bin_map[i] = []
|
||||
for i in builtins.range(len(bin_numbers)):
|
||||
bin_map[bin_numbers[i]].append(values[vv, i])
|
||||
return bin_map
|
||||
|
||||
|
||||
def _bin_edges(sample, bins=None, range=None):
|
||||
""" Create edge arrays
|
||||
"""
|
||||
Dlen, Ndim = sample.shape
|
||||
|
||||
nbin = np.empty(Ndim, int) # Number of bins in each dimension
|
||||
edges = Ndim * [None] # Bin edges for each dim (will be 2D array)
|
||||
dedges = Ndim * [None] # Spacing between edges (will be 2D array)
|
||||
|
||||
# Select range for each dimension
|
||||
# Used only if number of bins is given.
|
||||
if range is None:
|
||||
smin = np.atleast_1d(np.array(sample.min(axis=0), float))
|
||||
smax = np.atleast_1d(np.array(sample.max(axis=0), float))
|
||||
else:
|
||||
if len(range) != Ndim:
|
||||
raise ValueError(
|
||||
f"range given for {len(range)} dimensions; {Ndim} required")
|
||||
smin = np.empty(Ndim)
|
||||
smax = np.empty(Ndim)
|
||||
for i in builtins.range(Ndim):
|
||||
if range[i][1] < range[i][0]:
|
||||
raise ValueError(
|
||||
"In {}range, start must be <= stop".format(
|
||||
f"dimension {i + 1} of " if Ndim > 1 else ""))
|
||||
smin[i], smax[i] = range[i]
|
||||
|
||||
# Make sure the bins have a finite width.
|
||||
for i in builtins.range(len(smin)):
|
||||
if smin[i] == smax[i]:
|
||||
smin[i] = smin[i] - .5
|
||||
smax[i] = smax[i] + .5
|
||||
|
||||
# Preserve sample floating point precision in bin edges
|
||||
edges_dtype = (sample.dtype if np.issubdtype(sample.dtype, np.floating)
|
||||
else float)
|
||||
|
||||
# Create edge arrays
|
||||
for i in builtins.range(Ndim):
|
||||
if np.isscalar(bins[i]):
|
||||
nbin[i] = bins[i] + 2 # +2 for outlier bins
|
||||
edges[i] = np.linspace(smin[i], smax[i], nbin[i] - 1,
|
||||
dtype=edges_dtype)
|
||||
else:
|
||||
edges[i] = np.asarray(bins[i], edges_dtype)
|
||||
nbin[i] = len(edges[i]) + 1 # +1 for outlier bins
|
||||
dedges[i] = np.diff(edges[i])
|
||||
|
||||
nbin = np.asarray(nbin)
|
||||
|
||||
return nbin, edges, dedges
|
||||
|
||||
|
||||
def _bin_numbers(sample, nbin, edges, dedges):
|
||||
"""Compute the bin number each sample falls into, in each dimension
|
||||
"""
|
||||
Dlen, Ndim = sample.shape
|
||||
|
||||
sampBin = [
|
||||
np.digitize(sample[:, i], edges[i])
|
||||
for i in range(Ndim)
|
||||
]
|
||||
|
||||
# Using `digitize`, values that fall on an edge are put in the right bin.
|
||||
# For the rightmost bin, we want values equal to the right
|
||||
# edge to be counted in the last bin, and not as an outlier.
|
||||
for i in range(Ndim):
|
||||
# Find the rounding precision
|
||||
dedges_min = dedges[i].min()
|
||||
if dedges_min == 0:
|
||||
raise ValueError('The smallest edge difference is numerically 0.')
|
||||
decimal = int(-np.log10(dedges_min)) + 6
|
||||
# Find which points are on the rightmost edge.
|
||||
on_edge = np.where((sample[:, i] >= edges[i][-1]) &
|
||||
(np.around(sample[:, i], decimal) ==
|
||||
np.around(edges[i][-1], decimal)))[0]
|
||||
# Shift these points one bin to the left.
|
||||
sampBin[i][on_edge] -= 1
|
||||
|
||||
# Compute the sample indices in the flattened statistic matrix.
|
||||
binnumbers = np.ravel_multi_index(sampBin, nbin)
|
||||
|
||||
return binnumbers
|
||||
@@ -0,0 +1,375 @@
|
||||
from math import sqrt
|
||||
import numpy as np
|
||||
from scipy._lib._util import _validate_int
|
||||
from scipy.optimize import brentq
|
||||
from scipy.special import ndtri
|
||||
from ._discrete_distns import binom
|
||||
from ._common import ConfidenceInterval
|
||||
|
||||
|
||||
class BinomTestResult:
|
||||
"""
|
||||
Result of `scipy.stats.binomtest`.
|
||||
|
||||
Attributes
|
||||
----------
|
||||
k : int
|
||||
The number of successes (copied from `binomtest` input).
|
||||
n : int
|
||||
The number of trials (copied from `binomtest` input).
|
||||
alternative : str
|
||||
Indicates the alternative hypothesis specified in the input
|
||||
to `binomtest`. It will be one of ``'two-sided'``, ``'greater'``,
|
||||
or ``'less'``.
|
||||
statistic: float
|
||||
The estimate of the proportion of successes.
|
||||
pvalue : float
|
||||
The p-value of the hypothesis test.
|
||||
|
||||
"""
|
||||
def __init__(self, k, n, alternative, statistic, pvalue):
|
||||
self.k = k
|
||||
self.n = n
|
||||
self.alternative = alternative
|
||||
self.statistic = statistic
|
||||
self.pvalue = pvalue
|
||||
|
||||
# add alias for backward compatibility
|
||||
self.proportion_estimate = statistic
|
||||
|
||||
def __repr__(self):
|
||||
s = ("BinomTestResult("
|
||||
f"k={self.k}, "
|
||||
f"n={self.n}, "
|
||||
f"alternative={self.alternative!r}, "
|
||||
f"statistic={self.statistic}, "
|
||||
f"pvalue={self.pvalue})")
|
||||
return s
|
||||
|
||||
def proportion_ci(self, confidence_level=0.95, method='exact'):
|
||||
"""
|
||||
Compute the confidence interval for ``statistic``.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
confidence_level : float, optional
|
||||
Confidence level for the computed confidence interval
|
||||
of the estimated proportion. Default is 0.95.
|
||||
method : {'exact', 'wilson', 'wilsoncc'}, optional
|
||||
Selects the method used to compute the confidence interval
|
||||
for the estimate of the proportion:
|
||||
|
||||
'exact' :
|
||||
Use the Clopper-Pearson exact method [1]_.
|
||||
'wilson' :
|
||||
Wilson's method, without continuity correction ([2]_, [3]_).
|
||||
'wilsoncc' :
|
||||
Wilson's method, with continuity correction ([2]_, [3]_).
|
||||
|
||||
Default is ``'exact'``.
|
||||
|
||||
Returns
|
||||
-------
|
||||
ci : ``ConfidenceInterval`` object
|
||||
The object has attributes ``low`` and ``high`` that hold the
|
||||
lower and upper bounds of the confidence interval.
|
||||
|
||||
References
|
||||
----------
|
||||
.. [1] C. J. Clopper and E. S. Pearson, The use of confidence or
|
||||
fiducial limits illustrated in the case of the binomial,
|
||||
Biometrika, Vol. 26, No. 4, pp 404-413 (Dec. 1934).
|
||||
.. [2] E. B. Wilson, Probable inference, the law of succession, and
|
||||
statistical inference, J. Amer. Stat. Assoc., 22, pp 209-212
|
||||
(1927).
|
||||
.. [3] Robert G. Newcombe, Two-sided confidence intervals for the
|
||||
single proportion: comparison of seven methods, Statistics
|
||||
in Medicine, 17, pp 857-872 (1998).
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> from scipy.stats import binomtest
|
||||
>>> result = binomtest(k=7, n=50, p=0.1)
|
||||
>>> result.statistic
|
||||
0.14
|
||||
>>> result.proportion_ci()
|
||||
ConfidenceInterval(low=0.05819170033997342, high=0.26739600249700846)
|
||||
"""
|
||||
if method not in ('exact', 'wilson', 'wilsoncc'):
|
||||
raise ValueError(f"method ('{method}') must be one of 'exact', "
|
||||
"'wilson' or 'wilsoncc'.")
|
||||
if not (0 <= confidence_level <= 1):
|
||||
raise ValueError(f'confidence_level ({confidence_level}) must be in '
|
||||
'the interval [0, 1].')
|
||||
if method == 'exact':
|
||||
low, high = _binom_exact_conf_int(self.k, self.n,
|
||||
confidence_level,
|
||||
self.alternative)
|
||||
else:
|
||||
# method is 'wilson' or 'wilsoncc'
|
||||
low, high = _binom_wilson_conf_int(self.k, self.n,
|
||||
confidence_level,
|
||||
self.alternative,
|
||||
correction=method == 'wilsoncc')
|
||||
return ConfidenceInterval(low=low, high=high)
|
||||
|
||||
|
||||
def _findp(func):
|
||||
try:
|
||||
p = brentq(func, 0, 1)
|
||||
except RuntimeError:
|
||||
raise RuntimeError('numerical solver failed to converge when '
|
||||
'computing the confidence limits') from None
|
||||
except ValueError as exc:
|
||||
raise ValueError('brentq raised a ValueError; report this to the '
|
||||
'SciPy developers') from exc
|
||||
return p
|
||||
|
||||
|
||||
def _binom_exact_conf_int(k, n, confidence_level, alternative):
|
||||
"""
|
||||
Compute the estimate and confidence interval for the binomial test.
|
||||
|
||||
Returns proportion, prop_low, prop_high
|
||||
"""
|
||||
if alternative == 'two-sided':
|
||||
alpha = (1 - confidence_level) / 2
|
||||
if k == 0:
|
||||
plow = 0.0
|
||||
else:
|
||||
plow = _findp(lambda p: binom.sf(k-1, n, p) - alpha)
|
||||
if k == n:
|
||||
phigh = 1.0
|
||||
else:
|
||||
phigh = _findp(lambda p: binom.cdf(k, n, p) - alpha)
|
||||
elif alternative == 'less':
|
||||
alpha = 1 - confidence_level
|
||||
plow = 0.0
|
||||
if k == n:
|
||||
phigh = 1.0
|
||||
else:
|
||||
phigh = _findp(lambda p: binom.cdf(k, n, p) - alpha)
|
||||
elif alternative == 'greater':
|
||||
alpha = 1 - confidence_level
|
||||
if k == 0:
|
||||
plow = 0.0
|
||||
else:
|
||||
plow = _findp(lambda p: binom.sf(k-1, n, p) - alpha)
|
||||
phigh = 1.0
|
||||
return plow, phigh
|
||||
|
||||
|
||||
def _binom_wilson_conf_int(k, n, confidence_level, alternative, correction):
|
||||
# This function assumes that the arguments have already been validated.
|
||||
# In particular, `alternative` must be one of 'two-sided', 'less' or
|
||||
# 'greater'.
|
||||
p = k / n
|
||||
if alternative == 'two-sided':
|
||||
z = ndtri(0.5 + 0.5*confidence_level)
|
||||
else:
|
||||
z = ndtri(confidence_level)
|
||||
|
||||
# For reference, the formulas implemented here are from
|
||||
# Newcombe (1998) (ref. [3] in the proportion_ci docstring).
|
||||
denom = 2*(n + z**2)
|
||||
center = (2*n*p + z**2)/denom
|
||||
q = 1 - p
|
||||
if correction:
|
||||
if alternative == 'less' or k == 0:
|
||||
lo = 0.0
|
||||
else:
|
||||
dlo = (1 + z*sqrt(z**2 - 2 - 1/n + 4*p*(n*q + 1))) / denom
|
||||
lo = center - dlo
|
||||
if alternative == 'greater' or k == n:
|
||||
hi = 1.0
|
||||
else:
|
||||
dhi = (1 + z*sqrt(z**2 + 2 - 1/n + 4*p*(n*q - 1))) / denom
|
||||
hi = center + dhi
|
||||
else:
|
||||
delta = z/denom * sqrt(4*n*p*q + z**2)
|
||||
if alternative == 'less' or k == 0:
|
||||
lo = 0.0
|
||||
else:
|
||||
lo = center - delta
|
||||
if alternative == 'greater' or k == n:
|
||||
hi = 1.0
|
||||
else:
|
||||
hi = center + delta
|
||||
|
||||
return lo, hi
|
||||
|
||||
|
||||
def binomtest(k, n, p=0.5, alternative='two-sided'):
|
||||
"""
|
||||
Perform a test that the probability of success is p.
|
||||
|
||||
The binomial test [1]_ is a test of the null hypothesis that the
|
||||
probability of success in a Bernoulli experiment is `p`.
|
||||
|
||||
Details of the test can be found in many texts on statistics, such
|
||||
as section 24.5 of [2]_.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
k : int
|
||||
The number of successes.
|
||||
n : int
|
||||
The number of trials.
|
||||
p : float, optional
|
||||
The hypothesized probability of success, i.e. the expected
|
||||
proportion of successes. The value must be in the interval
|
||||
``0 <= p <= 1``. The default value is ``p = 0.5``.
|
||||
alternative : {'two-sided', 'greater', 'less'}, optional
|
||||
Indicates the alternative hypothesis. The default value is
|
||||
'two-sided'.
|
||||
|
||||
Returns
|
||||
-------
|
||||
result : `~scipy.stats._result_classes.BinomTestResult` instance
|
||||
The return value is an object with the following attributes:
|
||||
|
||||
k : int
|
||||
The number of successes (copied from `binomtest` input).
|
||||
n : int
|
||||
The number of trials (copied from `binomtest` input).
|
||||
alternative : str
|
||||
Indicates the alternative hypothesis specified in the input
|
||||
to `binomtest`. It will be one of ``'two-sided'``, ``'greater'``,
|
||||
or ``'less'``.
|
||||
statistic : float
|
||||
The estimate of the proportion of successes.
|
||||
pvalue : float
|
||||
The p-value of the hypothesis test.
|
||||
|
||||
The object has the following methods:
|
||||
|
||||
proportion_ci(confidence_level=0.95, method='exact') :
|
||||
Compute the confidence interval for ``statistic``.
|
||||
|
||||
Notes
|
||||
-----
|
||||
.. versionadded:: 1.7.0
|
||||
|
||||
References
|
||||
----------
|
||||
.. [1] Binomial test, https://en.wikipedia.org/wiki/Binomial_test
|
||||
.. [2] Jerrold H. Zar, Biostatistical Analysis (fifth edition),
|
||||
Prentice Hall, Upper Saddle River, New Jersey USA (2010)
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> from scipy.stats import binomtest
|
||||
|
||||
A car manufacturer claims that no more than 10% of their cars are unsafe.
|
||||
15 cars are inspected for safety, 3 were found to be unsafe. Test the
|
||||
manufacturer's claim:
|
||||
|
||||
>>> result = binomtest(3, n=15, p=0.1, alternative='greater')
|
||||
>>> result.pvalue
|
||||
0.18406106910639114
|
||||
|
||||
The null hypothesis cannot be rejected at the 5% level of significance
|
||||
because the returned p-value is greater than the critical value of 5%.
|
||||
|
||||
The test statistic is equal to the estimated proportion, which is simply
|
||||
``3/15``:
|
||||
|
||||
>>> result.statistic
|
||||
0.2
|
||||
|
||||
We can use the `proportion_ci()` method of the result to compute the
|
||||
confidence interval of the estimate:
|
||||
|
||||
>>> result.proportion_ci(confidence_level=0.95)
|
||||
ConfidenceInterval(low=0.05684686759024681, high=1.0)
|
||||
|
||||
"""
|
||||
k = _validate_int(k, 'k', minimum=0)
|
||||
n = _validate_int(n, 'n', minimum=1)
|
||||
if k > n:
|
||||
raise ValueError(f'k ({k}) must not be greater than n ({n}).')
|
||||
|
||||
if not (0 <= p <= 1):
|
||||
raise ValueError(f"p ({p}) must be in range [0,1]")
|
||||
|
||||
if alternative not in ('two-sided', 'less', 'greater'):
|
||||
raise ValueError(f"alternative ('{alternative}') not recognized; \n"
|
||||
"must be 'two-sided', 'less' or 'greater'")
|
||||
if alternative == 'less':
|
||||
pval = binom.cdf(k, n, p)
|
||||
elif alternative == 'greater':
|
||||
pval = binom.sf(k-1, n, p)
|
||||
else:
|
||||
# alternative is 'two-sided'
|
||||
d = binom.pmf(k, n, p)
|
||||
rerr = 1 + 1e-7
|
||||
if k == p * n:
|
||||
# special case as shortcut, would also be handled by `else` below
|
||||
pval = 1.
|
||||
elif k < p * n:
|
||||
ix = _binary_search_for_binom_tst(lambda x1: -binom.pmf(x1, n, p),
|
||||
-d*rerr, np.ceil(p * n), n)
|
||||
# y is the number of terms between mode and n that are <= d*rerr.
|
||||
# ix gave us the first term where a(ix) <= d*rerr < a(ix-1)
|
||||
# if the first equality doesn't hold, y=n-ix. Otherwise, we
|
||||
# need to include ix as well as the equality holds. Note that
|
||||
# the equality will hold in very very rare situations due to rerr.
|
||||
y = n - ix + int(d*rerr == binom.pmf(ix, n, p))
|
||||
pval = binom.cdf(k, n, p) + binom.sf(n - y, n, p)
|
||||
else:
|
||||
ix = _binary_search_for_binom_tst(lambda x1: binom.pmf(x1, n, p),
|
||||
d*rerr, 0, np.floor(p * n))
|
||||
# y is the number of terms between 0 and mode that are <= d*rerr.
|
||||
# we need to add a 1 to account for the 0 index.
|
||||
# For comparing this with old behavior, see
|
||||
# tst_binary_srch_for_binom_tst method in test_morestats.
|
||||
y = ix + 1
|
||||
pval = binom.cdf(y-1, n, p) + binom.sf(k-1, n, p)
|
||||
|
||||
pval = min(1.0, pval)
|
||||
|
||||
result = BinomTestResult(k=k, n=n, alternative=alternative,
|
||||
statistic=k/n, pvalue=pval)
|
||||
return result
|
||||
|
||||
|
||||
def _binary_search_for_binom_tst(a, d, lo, hi):
|
||||
"""
|
||||
Conducts an implicit binary search on a function specified by `a`.
|
||||
|
||||
Meant to be used on the binomial PMF for the case of two-sided tests
|
||||
to obtain the value on the other side of the mode where the tail
|
||||
probability should be computed. The values on either side of
|
||||
the mode are always in order, meaning binary search is applicable.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
a : callable
|
||||
The function over which to perform binary search. Its values
|
||||
for inputs lo and hi should be in ascending order.
|
||||
d : float
|
||||
The value to search.
|
||||
lo : int
|
||||
The lower end of range to search.
|
||||
hi : int
|
||||
The higher end of the range to search.
|
||||
|
||||
Returns
|
||||
-------
|
||||
int
|
||||
The index, i between lo and hi
|
||||
such that a(i)<=d<a(i+1)
|
||||
"""
|
||||
while lo < hi:
|
||||
mid = lo + (hi-lo)//2
|
||||
midval = a(mid)
|
||||
if midval < d:
|
||||
lo = mid+1
|
||||
elif midval > d:
|
||||
hi = mid-1
|
||||
else:
|
||||
return mid
|
||||
if a(lo) <= d:
|
||||
return lo
|
||||
else:
|
||||
return lo-1
|
||||
@@ -0,0 +1,53 @@
|
||||
from scipy.stats._boost.beta_ufunc import (
|
||||
_beta_pdf, _beta_cdf, _beta_sf, _beta_ppf,
|
||||
_beta_isf, _beta_mean, _beta_variance,
|
||||
_beta_skewness, _beta_kurtosis_excess,
|
||||
)
|
||||
|
||||
from scipy.stats._boost.binom_ufunc import (
|
||||
_binom_pdf, _binom_cdf, _binom_sf, _binom_ppf,
|
||||
_binom_isf, _binom_mean, _binom_variance,
|
||||
_binom_skewness, _binom_kurtosis_excess,
|
||||
)
|
||||
|
||||
from scipy.stats._boost.nbinom_ufunc import (
|
||||
_nbinom_pdf, _nbinom_cdf, _nbinom_sf, _nbinom_ppf,
|
||||
_nbinom_isf, _nbinom_mean, _nbinom_variance,
|
||||
_nbinom_skewness, _nbinom_kurtosis_excess,
|
||||
)
|
||||
|
||||
from scipy.stats._boost.hypergeom_ufunc import (
|
||||
_hypergeom_pdf, _hypergeom_cdf, _hypergeom_sf, _hypergeom_ppf,
|
||||
_hypergeom_isf, _hypergeom_mean, _hypergeom_variance,
|
||||
_hypergeom_skewness, _hypergeom_kurtosis_excess,
|
||||
)
|
||||
|
||||
from scipy.stats._boost.ncf_ufunc import (
|
||||
_ncf_pdf, _ncf_cdf, _ncf_sf, _ncf_ppf,
|
||||
_ncf_isf, _ncf_mean, _ncf_variance,
|
||||
_ncf_skewness, _ncf_kurtosis_excess,
|
||||
)
|
||||
|
||||
from scipy.stats._boost.ncx2_ufunc import (
|
||||
_ncx2_pdf, _ncx2_cdf, _ncx2_sf, _ncx2_ppf,
|
||||
_ncx2_isf, _ncx2_mean, _ncx2_variance,
|
||||
_ncx2_skewness, _ncx2_kurtosis_excess,
|
||||
)
|
||||
|
||||
from scipy.stats._boost.nct_ufunc import (
|
||||
_nct_pdf, _nct_cdf, _nct_sf, _nct_ppf,
|
||||
_nct_isf, _nct_mean, _nct_variance,
|
||||
_nct_skewness, _nct_kurtosis_excess,
|
||||
)
|
||||
|
||||
from scipy.stats._boost.skewnorm_ufunc import (
|
||||
_skewnorm_pdf, _skewnorm_cdf, _skewnorm_sf, _skewnorm_ppf,
|
||||
_skewnorm_isf, _skewnorm_mean, _skewnorm_variance,
|
||||
_skewnorm_skewness, _skewnorm_kurtosis_excess,
|
||||
)
|
||||
|
||||
from scipy.stats._boost.invgauss_ufunc import (
|
||||
_invgauss_pdf, _invgauss_cdf, _invgauss_sf, _invgauss_ppf,
|
||||
_invgauss_isf, _invgauss_mean, _invgauss_variance,
|
||||
_invgauss_skewness, _invgauss_kurtosis_excess,
|
||||
)
|
||||
BIN
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
BIN
Binary file not shown.
Binary file not shown.
BIN
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
BIN
Binary file not shown.
Binary file not shown.
@@ -0,0 +1,177 @@
|
||||
import numpy as np
|
||||
from functools import partial
|
||||
from scipy import stats
|
||||
|
||||
|
||||
def _bws_input_validation(x, y, alternative, method):
|
||||
''' Input validation and standardization for bws test'''
|
||||
x, y = np.atleast_1d(x, y)
|
||||
if x.ndim > 1 or y.ndim > 1:
|
||||
raise ValueError('`x` and `y` must be exactly one-dimensional.')
|
||||
if np.isnan(x).any() or np.isnan(y).any():
|
||||
raise ValueError('`x` and `y` must not contain NaNs.')
|
||||
if np.size(x) == 0 or np.size(y) == 0:
|
||||
raise ValueError('`x` and `y` must be of nonzero size.')
|
||||
|
||||
z = stats.rankdata(np.concatenate((x, y)))
|
||||
x, y = z[:len(x)], z[len(x):]
|
||||
|
||||
alternatives = {'two-sided', 'less', 'greater'}
|
||||
alternative = alternative.lower()
|
||||
if alternative not in alternatives:
|
||||
raise ValueError(f'`alternative` must be one of {alternatives}.')
|
||||
|
||||
method = stats.PermutationMethod() if method is None else method
|
||||
if not isinstance(method, stats.PermutationMethod):
|
||||
raise ValueError('`method` must be an instance of '
|
||||
'`scipy.stats.PermutationMethod`')
|
||||
|
||||
return x, y, alternative, method
|
||||
|
||||
|
||||
def _bws_statistic(x, y, alternative, axis):
|
||||
'''Compute the BWS test statistic for two independent samples'''
|
||||
# Public function currently does not accept `axis`, but `permutation_test`
|
||||
# uses `axis` to make vectorized call.
|
||||
|
||||
Ri, Hj = np.sort(x, axis=axis), np.sort(y, axis=axis)
|
||||
n, m = Ri.shape[axis], Hj.shape[axis]
|
||||
i, j = np.arange(1, n+1), np.arange(1, m+1)
|
||||
|
||||
Bx_num = Ri - (m + n)/n * i
|
||||
By_num = Hj - (m + n)/m * j
|
||||
|
||||
if alternative == 'two-sided':
|
||||
Bx_num *= Bx_num
|
||||
By_num *= By_num
|
||||
else:
|
||||
Bx_num *= np.abs(Bx_num)
|
||||
By_num *= np.abs(By_num)
|
||||
|
||||
Bx_den = i/(n+1) * (1 - i/(n+1)) * m*(m+n)/n
|
||||
By_den = j/(m+1) * (1 - j/(m+1)) * n*(m+n)/m
|
||||
|
||||
Bx = 1/n * np.sum(Bx_num/Bx_den, axis=axis)
|
||||
By = 1/m * np.sum(By_num/By_den, axis=axis)
|
||||
|
||||
B = (Bx + By) / 2 if alternative == 'two-sided' else (Bx - By) / 2
|
||||
|
||||
return B
|
||||
|
||||
|
||||
def bws_test(x, y, *, alternative="two-sided", method=None):
|
||||
r'''Perform the Baumgartner-Weiss-Schindler test on two independent samples.
|
||||
|
||||
The Baumgartner-Weiss-Schindler (BWS) test is a nonparametric test of
|
||||
the null hypothesis that the distribution underlying sample `x`
|
||||
is the same as the distribution underlying sample `y`. Unlike
|
||||
the Kolmogorov-Smirnov, Wilcoxon, and Cramer-Von Mises tests,
|
||||
the BWS test weights the integral by the variance of the difference
|
||||
in cumulative distribution functions (CDFs), emphasizing the tails of the
|
||||
distributions, which increases the power of the test in many applications.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x, y : array-like
|
||||
1-d arrays of samples.
|
||||
alternative : {'two-sided', 'less', 'greater'}, optional
|
||||
Defines the alternative hypothesis. Default is 'two-sided'.
|
||||
Let *F(u)* and *G(u)* be the cumulative distribution functions of the
|
||||
distributions underlying `x` and `y`, respectively. Then the following
|
||||
alternative hypotheses are available:
|
||||
|
||||
* 'two-sided': the distributions are not equal, i.e. *F(u) ≠ G(u)* for
|
||||
at least one *u*.
|
||||
* 'less': the distribution underlying `x` is stochastically less than
|
||||
the distribution underlying `y`, i.e. *F(u) >= G(u)* for all *u*.
|
||||
* 'greater': the distribution underlying `x` is stochastically greater
|
||||
than the distribution underlying `y`, i.e. *F(u) <= G(u)* for all
|
||||
*u*.
|
||||
|
||||
Under a more restrictive set of assumptions, the alternative hypotheses
|
||||
can be expressed in terms of the locations of the distributions;
|
||||
see [2] section 5.1.
|
||||
method : PermutationMethod, optional
|
||||
Configures the method used to compute the p-value. The default is
|
||||
the default `PermutationMethod` object.
|
||||
|
||||
Returns
|
||||
-------
|
||||
res : PermutationTestResult
|
||||
An object with attributes:
|
||||
|
||||
statistic : float
|
||||
The observed test statistic of the data.
|
||||
pvalue : float
|
||||
The p-value for the given alternative.
|
||||
null_distribution : ndarray
|
||||
The values of the test statistic generated under the null hypothesis.
|
||||
|
||||
See also
|
||||
--------
|
||||
scipy.stats.wilcoxon, scipy.stats.mannwhitneyu, scipy.stats.ttest_ind
|
||||
|
||||
Notes
|
||||
-----
|
||||
When ``alternative=='two-sided'``, the statistic is defined by the
|
||||
equations given in [1]_ Section 2. This statistic is not appropriate for
|
||||
one-sided alternatives; in that case, the statistic is the *negative* of
|
||||
that given by the equations in [1]_ Section 2. Consequently, when the
|
||||
distribution of the first sample is stochastically greater than that of the
|
||||
second sample, the statistic will tend to be positive.
|
||||
|
||||
References
|
||||
----------
|
||||
.. [1] Neuhäuser, M. (2005). Exact Tests Based on the
|
||||
Baumgartner-Weiss-Schindler Statistic: A Survey. Statistical Papers,
|
||||
46(1), 1-29.
|
||||
.. [2] Fay, M. P., & Proschan, M. A. (2010). Wilcoxon-Mann-Whitney or t-test?
|
||||
On assumptions for hypothesis tests and multiple interpretations of
|
||||
decision rules. Statistics surveys, 4, 1.
|
||||
|
||||
Examples
|
||||
--------
|
||||
We follow the example of table 3 in [1]_: Fourteen children were divided
|
||||
randomly into two groups. Their ranks at performing a specific tests are
|
||||
as follows.
|
||||
|
||||
>>> import numpy as np
|
||||
>>> x = [1, 2, 3, 4, 6, 7, 8]
|
||||
>>> y = [5, 9, 10, 11, 12, 13, 14]
|
||||
|
||||
We use the BWS test to assess whether there is a statistically significant
|
||||
difference between the two groups.
|
||||
The null hypothesis is that there is no difference in the distributions of
|
||||
performance between the two groups. We decide that a significance level of
|
||||
1% is required to reject the null hypothesis in favor of the alternative
|
||||
that the distributions are different.
|
||||
Since the number of samples is very small, we can compare the observed test
|
||||
statistic against the *exact* distribution of the test statistic under the
|
||||
null hypothesis.
|
||||
|
||||
>>> from scipy.stats import bws_test
|
||||
>>> res = bws_test(x, y)
|
||||
>>> print(res.statistic)
|
||||
5.132167152575315
|
||||
|
||||
This agrees with :math:`B = 5.132` reported in [1]_. The *p*-value produced
|
||||
by `bws_test` also agrees with :math:`p = 0.0029` reported in [1]_.
|
||||
|
||||
>>> print(res.pvalue)
|
||||
0.002913752913752914
|
||||
|
||||
Because the p-value is below our threshold of 1%, we take this as evidence
|
||||
against the null hypothesis in favor of the alternative that there is a
|
||||
difference in performance between the two groups.
|
||||
'''
|
||||
|
||||
x, y, alternative, method = _bws_input_validation(x, y, alternative,
|
||||
method)
|
||||
bws_statistic = partial(_bws_statistic, alternative=alternative)
|
||||
|
||||
permutation_alternative = 'less' if alternative == 'less' else 'greater'
|
||||
res = stats.permutation_test((x, y), bws_statistic,
|
||||
alternative=permutation_alternative,
|
||||
**method._asdict())
|
||||
|
||||
return res
|
||||
@@ -0,0 +1,459 @@
|
||||
import numpy as np
|
||||
|
||||
|
||||
def _validate_1d(a, name, allow_inf=False):
|
||||
if np.ndim(a) != 1:
|
||||
raise ValueError(f'`{name}` must be a one-dimensional sequence.')
|
||||
if np.isnan(a).any():
|
||||
raise ValueError(f'`{name}` must not contain nan.')
|
||||
if not allow_inf and np.isinf(a).any():
|
||||
raise ValueError(f'`{name}` must contain only finite values.')
|
||||
|
||||
|
||||
def _validate_interval(interval):
|
||||
interval = np.asarray(interval)
|
||||
if interval.shape == (0,):
|
||||
# The input was a sequence with length 0.
|
||||
interval = interval.reshape((0, 2))
|
||||
if interval.ndim != 2 or interval.shape[-1] != 2:
|
||||
raise ValueError('`interval` must be a two-dimensional array with '
|
||||
'shape (m, 2), where m is the number of '
|
||||
'interval-censored values, but got shape '
|
||||
f'{interval.shape}')
|
||||
|
||||
if np.isnan(interval).any():
|
||||
raise ValueError('`interval` must not contain nan.')
|
||||
if np.isinf(interval).all(axis=1).any():
|
||||
raise ValueError('In each row in `interval`, both values must not'
|
||||
' be infinite.')
|
||||
if (interval[:, 0] > interval[:, 1]).any():
|
||||
raise ValueError('In each row of `interval`, the left value must not'
|
||||
' exceed the right value.')
|
||||
|
||||
uncensored_mask = interval[:, 0] == interval[:, 1]
|
||||
left_mask = np.isinf(interval[:, 0])
|
||||
right_mask = np.isinf(interval[:, 1])
|
||||
interval_mask = np.isfinite(interval).all(axis=1) & ~uncensored_mask
|
||||
|
||||
uncensored2 = interval[uncensored_mask, 0]
|
||||
left2 = interval[left_mask, 1]
|
||||
right2 = interval[right_mask, 0]
|
||||
interval2 = interval[interval_mask]
|
||||
|
||||
return uncensored2, left2, right2, interval2
|
||||
|
||||
|
||||
def _validate_x_censored(x, censored):
|
||||
x = np.asarray(x)
|
||||
if x.ndim != 1:
|
||||
raise ValueError('`x` must be one-dimensional.')
|
||||
censored = np.asarray(censored)
|
||||
if censored.ndim != 1:
|
||||
raise ValueError('`censored` must be one-dimensional.')
|
||||
if (~np.isfinite(x)).any():
|
||||
raise ValueError('`x` must not contain nan or inf.')
|
||||
if censored.size != x.size:
|
||||
raise ValueError('`x` and `censored` must have the same length.')
|
||||
return x, censored.astype(bool)
|
||||
|
||||
|
||||
class CensoredData:
|
||||
"""
|
||||
Instances of this class represent censored data.
|
||||
|
||||
Instances may be passed to the ``fit`` method of continuous
|
||||
univariate SciPy distributions for maximum likelihood estimation.
|
||||
The *only* method of the univariate continuous distributions that
|
||||
understands `CensoredData` is the ``fit`` method. An instance of
|
||||
`CensoredData` can not be passed to methods such as ``pdf`` and
|
||||
``cdf``.
|
||||
|
||||
An observation is said to be *censored* when the precise value is unknown,
|
||||
but it has a known upper and/or lower bound. The conventional terminology
|
||||
is:
|
||||
|
||||
* left-censored: an observation is below a certain value but it is
|
||||
unknown by how much.
|
||||
* right-censored: an observation is above a certain value but it is
|
||||
unknown by how much.
|
||||
* interval-censored: an observation lies somewhere on an interval between
|
||||
two values.
|
||||
|
||||
Left-, right-, and interval-censored data can be represented by
|
||||
`CensoredData`.
|
||||
|
||||
For convenience, the class methods ``left_censored`` and
|
||||
``right_censored`` are provided to create a `CensoredData`
|
||||
instance from a single one-dimensional array of measurements
|
||||
and a corresponding boolean array to indicate which measurements
|
||||
are censored. The class method ``interval_censored`` accepts two
|
||||
one-dimensional arrays that hold the lower and upper bounds of the
|
||||
intervals.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
uncensored : array_like, 1D
|
||||
Uncensored observations.
|
||||
left : array_like, 1D
|
||||
Left-censored observations.
|
||||
right : array_like, 1D
|
||||
Right-censored observations.
|
||||
interval : array_like, 2D, with shape (m, 2)
|
||||
Interval-censored observations. Each row ``interval[k, :]``
|
||||
represents the interval for the kth interval-censored observation.
|
||||
|
||||
Notes
|
||||
-----
|
||||
In the input array `interval`, the lower bound of the interval may
|
||||
be ``-inf``, and the upper bound may be ``inf``, but at least one must be
|
||||
finite. When the lower bound is ``-inf``, the row represents a left-
|
||||
censored observation, and when the upper bound is ``inf``, the row
|
||||
represents a right-censored observation. If the length of an interval
|
||||
is 0 (i.e. ``interval[k, 0] == interval[k, 1]``, the observation is
|
||||
treated as uncensored. So one can represent all the types of censored
|
||||
and uncensored data in ``interval``, but it is generally more convenient
|
||||
to use `uncensored`, `left` and `right` for uncensored, left-censored and
|
||||
right-censored observations, respectively.
|
||||
|
||||
Examples
|
||||
--------
|
||||
In the most general case, a censored data set may contain values that
|
||||
are left-censored, right-censored, interval-censored, and uncensored.
|
||||
For example, here we create a data set with five observations. Two
|
||||
are uncensored (values 1 and 1.5), one is a left-censored observation
|
||||
of 0, one is a right-censored observation of 10 and one is
|
||||
interval-censored in the interval [2, 3].
|
||||
|
||||
>>> import numpy as np
|
||||
>>> from scipy.stats import CensoredData
|
||||
>>> data = CensoredData(uncensored=[1, 1.5], left=[0], right=[10],
|
||||
... interval=[[2, 3]])
|
||||
>>> print(data)
|
||||
CensoredData(5 values: 2 not censored, 1 left-censored,
|
||||
1 right-censored, 1 interval-censored)
|
||||
|
||||
Equivalently,
|
||||
|
||||
>>> data = CensoredData(interval=[[1, 1],
|
||||
... [1.5, 1.5],
|
||||
... [-np.inf, 0],
|
||||
... [10, np.inf],
|
||||
... [2, 3]])
|
||||
>>> print(data)
|
||||
CensoredData(5 values: 2 not censored, 1 left-censored,
|
||||
1 right-censored, 1 interval-censored)
|
||||
|
||||
A common case is to have a mix of uncensored observations and censored
|
||||
observations that are all right-censored (or all left-censored). For
|
||||
example, consider an experiment in which six devices are started at
|
||||
various times and left running until they fail. Assume that time is
|
||||
measured in hours, and the experiment is stopped after 30 hours, even
|
||||
if all the devices have not failed by that time. We might end up with
|
||||
data such as this::
|
||||
|
||||
Device Start-time Fail-time Time-to-failure
|
||||
1 0 13 13
|
||||
2 2 24 22
|
||||
3 5 22 17
|
||||
4 8 23 15
|
||||
5 10 *** >20
|
||||
6 12 *** >18
|
||||
|
||||
Two of the devices had not failed when the experiment was stopped;
|
||||
the observations of the time-to-failure for these two devices are
|
||||
right-censored. We can represent this data with
|
||||
|
||||
>>> data = CensoredData(uncensored=[13, 22, 17, 15], right=[20, 18])
|
||||
>>> print(data)
|
||||
CensoredData(6 values: 4 not censored, 2 right-censored)
|
||||
|
||||
Alternatively, we can use the method `CensoredData.right_censored` to
|
||||
create a representation of this data. The time-to-failure observations
|
||||
are put the list ``ttf``. The ``censored`` list indicates which values
|
||||
in ``ttf`` are censored.
|
||||
|
||||
>>> ttf = [13, 22, 17, 15, 20, 18]
|
||||
>>> censored = [False, False, False, False, True, True]
|
||||
|
||||
Pass these lists to `CensoredData.right_censored` to create an
|
||||
instance of `CensoredData`.
|
||||
|
||||
>>> data = CensoredData.right_censored(ttf, censored)
|
||||
>>> print(data)
|
||||
CensoredData(6 values: 4 not censored, 2 right-censored)
|
||||
|
||||
If the input data is interval censored and already stored in two
|
||||
arrays, one holding the low end of the intervals and another
|
||||
holding the high ends, the class method ``interval_censored`` can
|
||||
be used to create the `CensoredData` instance.
|
||||
|
||||
This example creates an instance with four interval-censored values.
|
||||
The intervals are [10, 11], [0.5, 1], [2, 3], and [12.5, 13.5].
|
||||
|
||||
>>> a = [10, 0.5, 2, 12.5] # Low ends of the intervals
|
||||
>>> b = [11, 1.0, 3, 13.5] # High ends of the intervals
|
||||
>>> data = CensoredData.interval_censored(low=a, high=b)
|
||||
>>> print(data)
|
||||
CensoredData(4 values: 0 not censored, 4 interval-censored)
|
||||
|
||||
Finally, we create and censor some data from the `weibull_min`
|
||||
distribution, and then fit `weibull_min` to that data. We'll assume
|
||||
that the location parameter is known to be 0.
|
||||
|
||||
>>> from scipy.stats import weibull_min
|
||||
>>> rng = np.random.default_rng()
|
||||
|
||||
Create the random data set.
|
||||
|
||||
>>> x = weibull_min.rvs(2.5, loc=0, scale=30, size=250, random_state=rng)
|
||||
>>> x[x > 40] = 40 # Right-censor values greater or equal to 40.
|
||||
|
||||
Create the `CensoredData` instance with the `right_censored` method.
|
||||
The censored values are those where the value is 40.
|
||||
|
||||
>>> data = CensoredData.right_censored(x, x == 40)
|
||||
>>> print(data)
|
||||
CensoredData(250 values: 215 not censored, 35 right-censored)
|
||||
|
||||
35 values have been right-censored.
|
||||
|
||||
Fit `weibull_min` to the censored data. We expect to shape and scale
|
||||
to be approximately 2.5 and 30, respectively.
|
||||
|
||||
>>> weibull_min.fit(data, floc=0)
|
||||
(2.3575922823897315, 0, 30.40650074451254)
|
||||
|
||||
"""
|
||||
|
||||
def __init__(self, uncensored=None, *, left=None, right=None,
|
||||
interval=None):
|
||||
if uncensored is None:
|
||||
uncensored = []
|
||||
if left is None:
|
||||
left = []
|
||||
if right is None:
|
||||
right = []
|
||||
if interval is None:
|
||||
interval = np.empty((0, 2))
|
||||
|
||||
_validate_1d(uncensored, 'uncensored')
|
||||
_validate_1d(left, 'left')
|
||||
_validate_1d(right, 'right')
|
||||
uncensored2, left2, right2, interval2 = _validate_interval(interval)
|
||||
|
||||
self._uncensored = np.concatenate((uncensored, uncensored2))
|
||||
self._left = np.concatenate((left, left2))
|
||||
self._right = np.concatenate((right, right2))
|
||||
# Note that by construction, the private attribute _interval
|
||||
# will be a 2D array that contains only finite values representing
|
||||
# intervals with nonzero but finite length.
|
||||
self._interval = interval2
|
||||
|
||||
def __repr__(self):
|
||||
uncensored_str = " ".join(np.array_repr(self._uncensored).split())
|
||||
left_str = " ".join(np.array_repr(self._left).split())
|
||||
right_str = " ".join(np.array_repr(self._right).split())
|
||||
interval_str = " ".join(np.array_repr(self._interval).split())
|
||||
return (f"CensoredData(uncensored={uncensored_str}, left={left_str}, "
|
||||
f"right={right_str}, interval={interval_str})")
|
||||
|
||||
def __str__(self):
|
||||
num_nc = len(self._uncensored)
|
||||
num_lc = len(self._left)
|
||||
num_rc = len(self._right)
|
||||
num_ic = len(self._interval)
|
||||
n = num_nc + num_lc + num_rc + num_ic
|
||||
parts = [f'{num_nc} not censored']
|
||||
if num_lc > 0:
|
||||
parts.append(f'{num_lc} left-censored')
|
||||
if num_rc > 0:
|
||||
parts.append(f'{num_rc} right-censored')
|
||||
if num_ic > 0:
|
||||
parts.append(f'{num_ic} interval-censored')
|
||||
return f'CensoredData({n} values: ' + ', '.join(parts) + ')'
|
||||
|
||||
# This is not a complete implementation of the arithmetic operators.
|
||||
# All we need is subtracting a scalar and dividing by a scalar.
|
||||
|
||||
def __sub__(self, other):
|
||||
return CensoredData(uncensored=self._uncensored - other,
|
||||
left=self._left - other,
|
||||
right=self._right - other,
|
||||
interval=self._interval - other)
|
||||
|
||||
def __truediv__(self, other):
|
||||
return CensoredData(uncensored=self._uncensored / other,
|
||||
left=self._left / other,
|
||||
right=self._right / other,
|
||||
interval=self._interval / other)
|
||||
|
||||
def __len__(self):
|
||||
"""
|
||||
The number of values (censored and not censored).
|
||||
"""
|
||||
return (len(self._uncensored) + len(self._left) + len(self._right)
|
||||
+ len(self._interval))
|
||||
|
||||
def num_censored(self):
|
||||
"""
|
||||
Number of censored values.
|
||||
"""
|
||||
return len(self._left) + len(self._right) + len(self._interval)
|
||||
|
||||
@classmethod
|
||||
def right_censored(cls, x, censored):
|
||||
"""
|
||||
Create a `CensoredData` instance of right-censored data.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x : array_like
|
||||
`x` is the array of observed data or measurements.
|
||||
`x` must be a one-dimensional sequence of finite numbers.
|
||||
censored : array_like of bool
|
||||
`censored` must be a one-dimensional sequence of boolean
|
||||
values. If ``censored[k]`` is True, the corresponding value
|
||||
in `x` is right-censored. That is, the value ``x[k]``
|
||||
is the lower bound of the true (but unknown) value.
|
||||
|
||||
Returns
|
||||
-------
|
||||
data : `CensoredData`
|
||||
An instance of `CensoredData` that represents the
|
||||
collection of uncensored and right-censored values.
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> from scipy.stats import CensoredData
|
||||
|
||||
Two uncensored values (4 and 10) and two right-censored values
|
||||
(24 and 25).
|
||||
|
||||
>>> data = CensoredData.right_censored([4, 10, 24, 25],
|
||||
... [False, False, True, True])
|
||||
>>> data
|
||||
CensoredData(uncensored=array([ 4., 10.]),
|
||||
left=array([], dtype=float64), right=array([24., 25.]),
|
||||
interval=array([], shape=(0, 2), dtype=float64))
|
||||
>>> print(data)
|
||||
CensoredData(4 values: 2 not censored, 2 right-censored)
|
||||
"""
|
||||
x, censored = _validate_x_censored(x, censored)
|
||||
return cls(uncensored=x[~censored], right=x[censored])
|
||||
|
||||
@classmethod
|
||||
def left_censored(cls, x, censored):
|
||||
"""
|
||||
Create a `CensoredData` instance of left-censored data.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x : array_like
|
||||
`x` is the array of observed data or measurements.
|
||||
`x` must be a one-dimensional sequence of finite numbers.
|
||||
censored : array_like of bool
|
||||
`censored` must be a one-dimensional sequence of boolean
|
||||
values. If ``censored[k]`` is True, the corresponding value
|
||||
in `x` is left-censored. That is, the value ``x[k]``
|
||||
is the upper bound of the true (but unknown) value.
|
||||
|
||||
Returns
|
||||
-------
|
||||
data : `CensoredData`
|
||||
An instance of `CensoredData` that represents the
|
||||
collection of uncensored and left-censored values.
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> from scipy.stats import CensoredData
|
||||
|
||||
Two uncensored values (0.12 and 0.033) and two left-censored values
|
||||
(both 1e-3).
|
||||
|
||||
>>> data = CensoredData.left_censored([0.12, 0.033, 1e-3, 1e-3],
|
||||
... [False, False, True, True])
|
||||
>>> data
|
||||
CensoredData(uncensored=array([0.12 , 0.033]),
|
||||
left=array([0.001, 0.001]), right=array([], dtype=float64),
|
||||
interval=array([], shape=(0, 2), dtype=float64))
|
||||
>>> print(data)
|
||||
CensoredData(4 values: 2 not censored, 2 left-censored)
|
||||
"""
|
||||
x, censored = _validate_x_censored(x, censored)
|
||||
return cls(uncensored=x[~censored], left=x[censored])
|
||||
|
||||
@classmethod
|
||||
def interval_censored(cls, low, high):
|
||||
"""
|
||||
Create a `CensoredData` instance of interval-censored data.
|
||||
|
||||
This method is useful when all the data is interval-censored, and
|
||||
the low and high ends of the intervals are already stored in
|
||||
separate one-dimensional arrays.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
low : array_like
|
||||
The one-dimensional array containing the low ends of the
|
||||
intervals.
|
||||
high : array_like
|
||||
The one-dimensional array containing the high ends of the
|
||||
intervals.
|
||||
|
||||
Returns
|
||||
-------
|
||||
data : `CensoredData`
|
||||
An instance of `CensoredData` that represents the
|
||||
collection of censored values.
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> import numpy as np
|
||||
>>> from scipy.stats import CensoredData
|
||||
|
||||
``a`` and ``b`` are the low and high ends of a collection of
|
||||
interval-censored values.
|
||||
|
||||
>>> a = [0.5, 2.0, 3.0, 5.5]
|
||||
>>> b = [1.0, 2.5, 3.5, 7.0]
|
||||
>>> data = CensoredData.interval_censored(low=a, high=b)
|
||||
>>> print(data)
|
||||
CensoredData(4 values: 0 not censored, 4 interval-censored)
|
||||
"""
|
||||
_validate_1d(low, 'low', allow_inf=True)
|
||||
_validate_1d(high, 'high', allow_inf=True)
|
||||
if len(low) != len(high):
|
||||
raise ValueError('`low` and `high` must have the same length.')
|
||||
interval = np.column_stack((low, high))
|
||||
uncensored, left, right, interval = _validate_interval(interval)
|
||||
return cls(uncensored=uncensored, left=left, right=right,
|
||||
interval=interval)
|
||||
|
||||
def _uncensor(self):
|
||||
"""
|
||||
This function is used when a non-censored version of the data
|
||||
is needed to create a rough estimate of the parameters of a
|
||||
distribution via the method of moments or some similar method.
|
||||
The data is "uncensored" by taking the given endpoints as the
|
||||
data for the left- or right-censored data, and the mean for the
|
||||
interval-censored data.
|
||||
"""
|
||||
data = np.concatenate((self._uncensored, self._left, self._right,
|
||||
self._interval.mean(axis=1)))
|
||||
return data
|
||||
|
||||
def _supported(self, a, b):
|
||||
"""
|
||||
Return a subset of self containing the values that are in
|
||||
(or overlap with) the interval (a, b).
|
||||
"""
|
||||
uncensored = self._uncensored
|
||||
uncensored = uncensored[(a < uncensored) & (uncensored < b)]
|
||||
left = self._left
|
||||
left = left[a < left]
|
||||
right = self._right
|
||||
right = right[right < b]
|
||||
interval = self._interval
|
||||
interval = interval[(a < interval[:, 1]) & (interval[:, 0] < b)]
|
||||
return CensoredData(uncensored, left=left, right=right,
|
||||
interval=interval)
|
||||
@@ -0,0 +1,5 @@
|
||||
from collections import namedtuple
|
||||
|
||||
|
||||
ConfidenceInterval = namedtuple("ConfidenceInterval", ["low", "high"])
|
||||
ConfidenceInterval. __doc__ = "Class for confidence intervals."
|
||||
@@ -0,0 +1,39 @@
|
||||
"""
|
||||
Statistics-related constants.
|
||||
|
||||
"""
|
||||
import numpy as np
|
||||
|
||||
|
||||
# The smallest representable positive number such that 1.0 + _EPS != 1.0.
|
||||
_EPS = np.finfo(float).eps
|
||||
|
||||
# The largest [in magnitude] usable floating value.
|
||||
_XMAX = np.finfo(float).max
|
||||
|
||||
# The log of the largest usable floating value; useful for knowing
|
||||
# when exp(something) will overflow
|
||||
_LOGXMAX = np.log(_XMAX)
|
||||
|
||||
# The smallest [in magnitude] usable (i.e. not subnormal) double precision
|
||||
# floating value.
|
||||
_XMIN = np.finfo(float).tiny
|
||||
|
||||
# The log of the smallest [in magnitude] usable (i.e not subnormal)
|
||||
# double precision floating value.
|
||||
_LOGXMIN = np.log(_XMIN)
|
||||
|
||||
# -special.psi(1)
|
||||
_EULER = 0.577215664901532860606512090082402431042
|
||||
|
||||
# special.zeta(3, 1) Apery's constant
|
||||
_ZETA3 = 1.202056903159594285399738161511449990765
|
||||
|
||||
# sqrt(pi)
|
||||
_SQRT_PI = 1.772453850905516027298167483341145182798
|
||||
|
||||
# sqrt(2/pi)
|
||||
_SQRT_2_OVER_PI = 0.7978845608028654
|
||||
|
||||
# log(sqrt(2/pi))
|
||||
_LOG_SQRT_2_OVER_PI = -0.22579135264472744
|
||||
File diff suppressed because it is too large
Load Diff
@@ -0,0 +1,633 @@
|
||||
from functools import cached_property
|
||||
|
||||
import numpy as np
|
||||
from scipy import linalg
|
||||
from scipy.stats import _multivariate
|
||||
|
||||
|
||||
__all__ = ["Covariance"]
|
||||
|
||||
|
||||
class Covariance:
|
||||
"""
|
||||
Representation of a covariance matrix
|
||||
|
||||
Calculations involving covariance matrices (e.g. data whitening,
|
||||
multivariate normal function evaluation) are often performed more
|
||||
efficiently using a decomposition of the covariance matrix instead of the
|
||||
covariance matrix itself. This class allows the user to construct an
|
||||
object representing a covariance matrix using any of several
|
||||
decompositions and perform calculations using a common interface.
|
||||
|
||||
.. note::
|
||||
|
||||
The `Covariance` class cannot be instantiated directly. Instead, use
|
||||
one of the factory methods (e.g. `Covariance.from_diagonal`).
|
||||
|
||||
Examples
|
||||
--------
|
||||
The `Covariance` class is is used by calling one of its
|
||||
factory methods to create a `Covariance` object, then pass that
|
||||
representation of the `Covariance` matrix as a shape parameter of a
|
||||
multivariate distribution.
|
||||
|
||||
For instance, the multivariate normal distribution can accept an array
|
||||
representing a covariance matrix:
|
||||
|
||||
>>> from scipy import stats
|
||||
>>> import numpy as np
|
||||
>>> d = [1, 2, 3]
|
||||
>>> A = np.diag(d) # a diagonal covariance matrix
|
||||
>>> x = [4, -2, 5] # a point of interest
|
||||
>>> dist = stats.multivariate_normal(mean=[0, 0, 0], cov=A)
|
||||
>>> dist.pdf(x)
|
||||
4.9595685102808205e-08
|
||||
|
||||
but the calculations are performed in a very generic way that does not
|
||||
take advantage of any special properties of the covariance matrix. Because
|
||||
our covariance matrix is diagonal, we can use ``Covariance.from_diagonal``
|
||||
to create an object representing the covariance matrix, and
|
||||
`multivariate_normal` can use this to compute the probability density
|
||||
function more efficiently.
|
||||
|
||||
>>> cov = stats.Covariance.from_diagonal(d)
|
||||
>>> dist = stats.multivariate_normal(mean=[0, 0, 0], cov=cov)
|
||||
>>> dist.pdf(x)
|
||||
4.9595685102808205e-08
|
||||
|
||||
"""
|
||||
def __init__(self):
|
||||
message = ("The `Covariance` class cannot be instantiated directly. "
|
||||
"Please use one of the factory methods "
|
||||
"(e.g. `Covariance.from_diagonal`).")
|
||||
raise NotImplementedError(message)
|
||||
|
||||
@staticmethod
|
||||
def from_diagonal(diagonal):
|
||||
r"""
|
||||
Return a representation of a covariance matrix from its diagonal.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
diagonal : array_like
|
||||
The diagonal elements of a diagonal matrix.
|
||||
|
||||
Notes
|
||||
-----
|
||||
Let the diagonal elements of a diagonal covariance matrix :math:`D` be
|
||||
stored in the vector :math:`d`.
|
||||
|
||||
When all elements of :math:`d` are strictly positive, whitening of a
|
||||
data point :math:`x` is performed by computing
|
||||
:math:`x \cdot d^{-1/2}`, where the inverse square root can be taken
|
||||
element-wise.
|
||||
:math:`\log\det{D}` is calculated as :math:`-2 \sum(\log{d})`,
|
||||
where the :math:`\log` operation is performed element-wise.
|
||||
|
||||
This `Covariance` class supports singular covariance matrices. When
|
||||
computing ``_log_pdet``, non-positive elements of :math:`d` are
|
||||
ignored. Whitening is not well defined when the point to be whitened
|
||||
does not lie in the span of the columns of the covariance matrix. The
|
||||
convention taken here is to treat the inverse square root of
|
||||
non-positive elements of :math:`d` as zeros.
|
||||
|
||||
Examples
|
||||
--------
|
||||
Prepare a symmetric positive definite covariance matrix ``A`` and a
|
||||
data point ``x``.
|
||||
|
||||
>>> import numpy as np
|
||||
>>> from scipy import stats
|
||||
>>> rng = np.random.default_rng()
|
||||
>>> n = 5
|
||||
>>> A = np.diag(rng.random(n))
|
||||
>>> x = rng.random(size=n)
|
||||
|
||||
Extract the diagonal from ``A`` and create the `Covariance` object.
|
||||
|
||||
>>> d = np.diag(A)
|
||||
>>> cov = stats.Covariance.from_diagonal(d)
|
||||
|
||||
Compare the functionality of the `Covariance` object against a
|
||||
reference implementations.
|
||||
|
||||
>>> res = cov.whiten(x)
|
||||
>>> ref = np.diag(d**-0.5) @ x
|
||||
>>> np.allclose(res, ref)
|
||||
True
|
||||
>>> res = cov.log_pdet
|
||||
>>> ref = np.linalg.slogdet(A)[-1]
|
||||
>>> np.allclose(res, ref)
|
||||
True
|
||||
|
||||
"""
|
||||
return CovViaDiagonal(diagonal)
|
||||
|
||||
@staticmethod
|
||||
def from_precision(precision, covariance=None):
|
||||
r"""
|
||||
Return a representation of a covariance from its precision matrix.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
precision : array_like
|
||||
The precision matrix; that is, the inverse of a square, symmetric,
|
||||
positive definite covariance matrix.
|
||||
covariance : array_like, optional
|
||||
The square, symmetric, positive definite covariance matrix. If not
|
||||
provided, this may need to be calculated (e.g. to evaluate the
|
||||
cumulative distribution function of
|
||||
`scipy.stats.multivariate_normal`) by inverting `precision`.
|
||||
|
||||
Notes
|
||||
-----
|
||||
Let the covariance matrix be :math:`A`, its precision matrix be
|
||||
:math:`P = A^{-1}`, and :math:`L` be the lower Cholesky factor such
|
||||
that :math:`L L^T = P`.
|
||||
Whitening of a data point :math:`x` is performed by computing
|
||||
:math:`x^T L`. :math:`\log\det{A}` is calculated as
|
||||
:math:`-2tr(\log{L})`, where the :math:`\log` operation is performed
|
||||
element-wise.
|
||||
|
||||
This `Covariance` class does not support singular covariance matrices
|
||||
because the precision matrix does not exist for a singular covariance
|
||||
matrix.
|
||||
|
||||
Examples
|
||||
--------
|
||||
Prepare a symmetric positive definite precision matrix ``P`` and a
|
||||
data point ``x``. (If the precision matrix is not already available,
|
||||
consider the other factory methods of the `Covariance` class.)
|
||||
|
||||
>>> import numpy as np
|
||||
>>> from scipy import stats
|
||||
>>> rng = np.random.default_rng()
|
||||
>>> n = 5
|
||||
>>> P = rng.random(size=(n, n))
|
||||
>>> P = P @ P.T # a precision matrix must be positive definite
|
||||
>>> x = rng.random(size=n)
|
||||
|
||||
Create the `Covariance` object.
|
||||
|
||||
>>> cov = stats.Covariance.from_precision(P)
|
||||
|
||||
Compare the functionality of the `Covariance` object against
|
||||
reference implementations.
|
||||
|
||||
>>> res = cov.whiten(x)
|
||||
>>> ref = x @ np.linalg.cholesky(P)
|
||||
>>> np.allclose(res, ref)
|
||||
True
|
||||
>>> res = cov.log_pdet
|
||||
>>> ref = -np.linalg.slogdet(P)[-1]
|
||||
>>> np.allclose(res, ref)
|
||||
True
|
||||
|
||||
"""
|
||||
return CovViaPrecision(precision, covariance)
|
||||
|
||||
@staticmethod
|
||||
def from_cholesky(cholesky):
|
||||
r"""
|
||||
Representation of a covariance provided via the (lower) Cholesky factor
|
||||
|
||||
Parameters
|
||||
----------
|
||||
cholesky : array_like
|
||||
The lower triangular Cholesky factor of the covariance matrix.
|
||||
|
||||
Notes
|
||||
-----
|
||||
Let the covariance matrix be :math:`A` and :math:`L` be the lower
|
||||
Cholesky factor such that :math:`L L^T = A`.
|
||||
Whitening of a data point :math:`x` is performed by computing
|
||||
:math:`L^{-1} x`. :math:`\log\det{A}` is calculated as
|
||||
:math:`2tr(\log{L})`, where the :math:`\log` operation is performed
|
||||
element-wise.
|
||||
|
||||
This `Covariance` class does not support singular covariance matrices
|
||||
because the Cholesky decomposition does not exist for a singular
|
||||
covariance matrix.
|
||||
|
||||
Examples
|
||||
--------
|
||||
Prepare a symmetric positive definite covariance matrix ``A`` and a
|
||||
data point ``x``.
|
||||
|
||||
>>> import numpy as np
|
||||
>>> from scipy import stats
|
||||
>>> rng = np.random.default_rng()
|
||||
>>> n = 5
|
||||
>>> A = rng.random(size=(n, n))
|
||||
>>> A = A @ A.T # make the covariance symmetric positive definite
|
||||
>>> x = rng.random(size=n)
|
||||
|
||||
Perform the Cholesky decomposition of ``A`` and create the
|
||||
`Covariance` object.
|
||||
|
||||
>>> L = np.linalg.cholesky(A)
|
||||
>>> cov = stats.Covariance.from_cholesky(L)
|
||||
|
||||
Compare the functionality of the `Covariance` object against
|
||||
reference implementation.
|
||||
|
||||
>>> from scipy.linalg import solve_triangular
|
||||
>>> res = cov.whiten(x)
|
||||
>>> ref = solve_triangular(L, x, lower=True)
|
||||
>>> np.allclose(res, ref)
|
||||
True
|
||||
>>> res = cov.log_pdet
|
||||
>>> ref = np.linalg.slogdet(A)[-1]
|
||||
>>> np.allclose(res, ref)
|
||||
True
|
||||
|
||||
"""
|
||||
return CovViaCholesky(cholesky)
|
||||
|
||||
@staticmethod
|
||||
def from_eigendecomposition(eigendecomposition):
|
||||
r"""
|
||||
Representation of a covariance provided via eigendecomposition
|
||||
|
||||
Parameters
|
||||
----------
|
||||
eigendecomposition : sequence
|
||||
A sequence (nominally a tuple) containing the eigenvalue and
|
||||
eigenvector arrays as computed by `scipy.linalg.eigh` or
|
||||
`numpy.linalg.eigh`.
|
||||
|
||||
Notes
|
||||
-----
|
||||
Let the covariance matrix be :math:`A`, let :math:`V` be matrix of
|
||||
eigenvectors, and let :math:`W` be the diagonal matrix of eigenvalues
|
||||
such that `V W V^T = A`.
|
||||
|
||||
When all of the eigenvalues are strictly positive, whitening of a
|
||||
data point :math:`x` is performed by computing
|
||||
:math:`x^T (V W^{-1/2})`, where the inverse square root can be taken
|
||||
element-wise.
|
||||
:math:`\log\det{A}` is calculated as :math:`tr(\log{W})`,
|
||||
where the :math:`\log` operation is performed element-wise.
|
||||
|
||||
This `Covariance` class supports singular covariance matrices. When
|
||||
computing ``_log_pdet``, non-positive eigenvalues are ignored.
|
||||
Whitening is not well defined when the point to be whitened
|
||||
does not lie in the span of the columns of the covariance matrix. The
|
||||
convention taken here is to treat the inverse square root of
|
||||
non-positive eigenvalues as zeros.
|
||||
|
||||
Examples
|
||||
--------
|
||||
Prepare a symmetric positive definite covariance matrix ``A`` and a
|
||||
data point ``x``.
|
||||
|
||||
>>> import numpy as np
|
||||
>>> from scipy import stats
|
||||
>>> rng = np.random.default_rng()
|
||||
>>> n = 5
|
||||
>>> A = rng.random(size=(n, n))
|
||||
>>> A = A @ A.T # make the covariance symmetric positive definite
|
||||
>>> x = rng.random(size=n)
|
||||
|
||||
Perform the eigendecomposition of ``A`` and create the `Covariance`
|
||||
object.
|
||||
|
||||
>>> w, v = np.linalg.eigh(A)
|
||||
>>> cov = stats.Covariance.from_eigendecomposition((w, v))
|
||||
|
||||
Compare the functionality of the `Covariance` object against
|
||||
reference implementations.
|
||||
|
||||
>>> res = cov.whiten(x)
|
||||
>>> ref = x @ (v @ np.diag(w**-0.5))
|
||||
>>> np.allclose(res, ref)
|
||||
True
|
||||
>>> res = cov.log_pdet
|
||||
>>> ref = np.linalg.slogdet(A)[-1]
|
||||
>>> np.allclose(res, ref)
|
||||
True
|
||||
|
||||
"""
|
||||
return CovViaEigendecomposition(eigendecomposition)
|
||||
|
||||
def whiten(self, x):
|
||||
"""
|
||||
Perform a whitening transformation on data.
|
||||
|
||||
"Whitening" ("white" as in "white noise", in which each frequency has
|
||||
equal magnitude) transforms a set of random variables into a new set of
|
||||
random variables with unit-diagonal covariance. When a whitening
|
||||
transform is applied to a sample of points distributed according to
|
||||
a multivariate normal distribution with zero mean, the covariance of
|
||||
the transformed sample is approximately the identity matrix.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x : array_like
|
||||
An array of points. The last dimension must correspond with the
|
||||
dimensionality of the space, i.e., the number of columns in the
|
||||
covariance matrix.
|
||||
|
||||
Returns
|
||||
-------
|
||||
x_ : array_like
|
||||
The transformed array of points.
|
||||
|
||||
References
|
||||
----------
|
||||
.. [1] "Whitening Transformation". Wikipedia.
|
||||
https://en.wikipedia.org/wiki/Whitening_transformation
|
||||
.. [2] Novak, Lukas, and Miroslav Vorechovsky. "Generalization of
|
||||
coloring linear transformation". Transactions of VSB 18.2
|
||||
(2018): 31-35. :doi:`10.31490/tces-2018-0013`
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> import numpy as np
|
||||
>>> from scipy import stats
|
||||
>>> rng = np.random.default_rng()
|
||||
>>> n = 3
|
||||
>>> A = rng.random(size=(n, n))
|
||||
>>> cov_array = A @ A.T # make matrix symmetric positive definite
|
||||
>>> precision = np.linalg.inv(cov_array)
|
||||
>>> cov_object = stats.Covariance.from_precision(precision)
|
||||
>>> x = rng.multivariate_normal(np.zeros(n), cov_array, size=(10000))
|
||||
>>> x_ = cov_object.whiten(x)
|
||||
>>> np.cov(x_, rowvar=False) # near-identity covariance
|
||||
array([[0.97862122, 0.00893147, 0.02430451],
|
||||
[0.00893147, 0.96719062, 0.02201312],
|
||||
[0.02430451, 0.02201312, 0.99206881]])
|
||||
|
||||
"""
|
||||
return self._whiten(np.asarray(x))
|
||||
|
||||
def colorize(self, x):
|
||||
"""
|
||||
Perform a colorizing transformation on data.
|
||||
|
||||
"Colorizing" ("color" as in "colored noise", in which different
|
||||
frequencies may have different magnitudes) transforms a set of
|
||||
uncorrelated random variables into a new set of random variables with
|
||||
the desired covariance. When a coloring transform is applied to a
|
||||
sample of points distributed according to a multivariate normal
|
||||
distribution with identity covariance and zero mean, the covariance of
|
||||
the transformed sample is approximately the covariance matrix used
|
||||
in the coloring transform.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x : array_like
|
||||
An array of points. The last dimension must correspond with the
|
||||
dimensionality of the space, i.e., the number of columns in the
|
||||
covariance matrix.
|
||||
|
||||
Returns
|
||||
-------
|
||||
x_ : array_like
|
||||
The transformed array of points.
|
||||
|
||||
References
|
||||
----------
|
||||
.. [1] "Whitening Transformation". Wikipedia.
|
||||
https://en.wikipedia.org/wiki/Whitening_transformation
|
||||
.. [2] Novak, Lukas, and Miroslav Vorechovsky. "Generalization of
|
||||
coloring linear transformation". Transactions of VSB 18.2
|
||||
(2018): 31-35. :doi:`10.31490/tces-2018-0013`
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> import numpy as np
|
||||
>>> from scipy import stats
|
||||
>>> rng = np.random.default_rng(1638083107694713882823079058616272161)
|
||||
>>> n = 3
|
||||
>>> A = rng.random(size=(n, n))
|
||||
>>> cov_array = A @ A.T # make matrix symmetric positive definite
|
||||
>>> cholesky = np.linalg.cholesky(cov_array)
|
||||
>>> cov_object = stats.Covariance.from_cholesky(cholesky)
|
||||
>>> x = rng.multivariate_normal(np.zeros(n), np.eye(n), size=(10000))
|
||||
>>> x_ = cov_object.colorize(x)
|
||||
>>> cov_data = np.cov(x_, rowvar=False)
|
||||
>>> np.allclose(cov_data, cov_array, rtol=3e-2)
|
||||
True
|
||||
"""
|
||||
return self._colorize(np.asarray(x))
|
||||
|
||||
@property
|
||||
def log_pdet(self):
|
||||
"""
|
||||
Log of the pseudo-determinant of the covariance matrix
|
||||
"""
|
||||
return np.array(self._log_pdet, dtype=float)[()]
|
||||
|
||||
@property
|
||||
def rank(self):
|
||||
"""
|
||||
Rank of the covariance matrix
|
||||
"""
|
||||
return np.array(self._rank, dtype=int)[()]
|
||||
|
||||
@property
|
||||
def covariance(self):
|
||||
"""
|
||||
Explicit representation of the covariance matrix
|
||||
"""
|
||||
return self._covariance
|
||||
|
||||
@property
|
||||
def shape(self):
|
||||
"""
|
||||
Shape of the covariance array
|
||||
"""
|
||||
return self._shape
|
||||
|
||||
def _validate_matrix(self, A, name):
|
||||
A = np.atleast_2d(A)
|
||||
m, n = A.shape[-2:]
|
||||
if m != n or A.ndim != 2 or not (np.issubdtype(A.dtype, np.integer) or
|
||||
np.issubdtype(A.dtype, np.floating)):
|
||||
message = (f"The input `{name}` must be a square, "
|
||||
"two-dimensional array of real numbers.")
|
||||
raise ValueError(message)
|
||||
return A
|
||||
|
||||
def _validate_vector(self, A, name):
|
||||
A = np.atleast_1d(A)
|
||||
if A.ndim != 1 or not (np.issubdtype(A.dtype, np.integer) or
|
||||
np.issubdtype(A.dtype, np.floating)):
|
||||
message = (f"The input `{name}` must be a one-dimensional array "
|
||||
"of real numbers.")
|
||||
raise ValueError(message)
|
||||
return A
|
||||
|
||||
|
||||
class CovViaPrecision(Covariance):
|
||||
|
||||
def __init__(self, precision, covariance=None):
|
||||
precision = self._validate_matrix(precision, 'precision')
|
||||
if covariance is not None:
|
||||
covariance = self._validate_matrix(covariance, 'covariance')
|
||||
message = "`precision.shape` must equal `covariance.shape`."
|
||||
if precision.shape != covariance.shape:
|
||||
raise ValueError(message)
|
||||
|
||||
self._chol_P = np.linalg.cholesky(precision)
|
||||
self._log_pdet = -2*np.log(np.diag(self._chol_P)).sum(axis=-1)
|
||||
self._rank = precision.shape[-1] # must be full rank if invertible
|
||||
self._precision = precision
|
||||
self._cov_matrix = covariance
|
||||
self._shape = precision.shape
|
||||
self._allow_singular = False
|
||||
|
||||
def _whiten(self, x):
|
||||
return x @ self._chol_P
|
||||
|
||||
@cached_property
|
||||
def _covariance(self):
|
||||
n = self._shape[-1]
|
||||
return (linalg.cho_solve((self._chol_P, True), np.eye(n))
|
||||
if self._cov_matrix is None else self._cov_matrix)
|
||||
|
||||
def _colorize(self, x):
|
||||
return linalg.solve_triangular(self._chol_P.T, x.T, lower=False).T
|
||||
|
||||
|
||||
def _dot_diag(x, d):
|
||||
# If d were a full diagonal matrix, x @ d would always do what we want.
|
||||
# Special treatment is needed for n-dimensional `d` in which each row
|
||||
# includes only the diagonal elements of a covariance matrix.
|
||||
return x * d if x.ndim < 2 else x * np.expand_dims(d, -2)
|
||||
|
||||
|
||||
class CovViaDiagonal(Covariance):
|
||||
|
||||
def __init__(self, diagonal):
|
||||
diagonal = self._validate_vector(diagonal, 'diagonal')
|
||||
|
||||
i_zero = diagonal <= 0
|
||||
positive_diagonal = np.array(diagonal, dtype=np.float64)
|
||||
|
||||
positive_diagonal[i_zero] = 1 # ones don't affect determinant
|
||||
self._log_pdet = np.sum(np.log(positive_diagonal), axis=-1)
|
||||
|
||||
psuedo_reciprocals = 1 / np.sqrt(positive_diagonal)
|
||||
psuedo_reciprocals[i_zero] = 0
|
||||
|
||||
self._sqrt_diagonal = np.sqrt(diagonal)
|
||||
self._LP = psuedo_reciprocals
|
||||
self._rank = positive_diagonal.shape[-1] - i_zero.sum(axis=-1)
|
||||
self._covariance = np.apply_along_axis(np.diag, -1, diagonal)
|
||||
self._i_zero = i_zero
|
||||
self._shape = self._covariance.shape
|
||||
self._allow_singular = True
|
||||
|
||||
def _whiten(self, x):
|
||||
return _dot_diag(x, self._LP)
|
||||
|
||||
def _colorize(self, x):
|
||||
return _dot_diag(x, self._sqrt_diagonal)
|
||||
|
||||
def _support_mask(self, x):
|
||||
"""
|
||||
Check whether x lies in the support of the distribution.
|
||||
"""
|
||||
return ~np.any(_dot_diag(x, self._i_zero), axis=-1)
|
||||
|
||||
|
||||
class CovViaCholesky(Covariance):
|
||||
|
||||
def __init__(self, cholesky):
|
||||
L = self._validate_matrix(cholesky, 'cholesky')
|
||||
|
||||
self._factor = L
|
||||
self._log_pdet = 2*np.log(np.diag(self._factor)).sum(axis=-1)
|
||||
self._rank = L.shape[-1] # must be full rank for cholesky
|
||||
self._shape = L.shape
|
||||
self._allow_singular = False
|
||||
|
||||
@cached_property
|
||||
def _covariance(self):
|
||||
return self._factor @ self._factor.T
|
||||
|
||||
def _whiten(self, x):
|
||||
res = linalg.solve_triangular(self._factor, x.T, lower=True).T
|
||||
return res
|
||||
|
||||
def _colorize(self, x):
|
||||
return x @ self._factor.T
|
||||
|
||||
|
||||
class CovViaEigendecomposition(Covariance):
|
||||
|
||||
def __init__(self, eigendecomposition):
|
||||
eigenvalues, eigenvectors = eigendecomposition
|
||||
eigenvalues = self._validate_vector(eigenvalues, 'eigenvalues')
|
||||
eigenvectors = self._validate_matrix(eigenvectors, 'eigenvectors')
|
||||
message = ("The shapes of `eigenvalues` and `eigenvectors` "
|
||||
"must be compatible.")
|
||||
try:
|
||||
eigenvalues = np.expand_dims(eigenvalues, -2)
|
||||
eigenvectors, eigenvalues = np.broadcast_arrays(eigenvectors,
|
||||
eigenvalues)
|
||||
eigenvalues = eigenvalues[..., 0, :]
|
||||
except ValueError:
|
||||
raise ValueError(message)
|
||||
|
||||
i_zero = eigenvalues <= 0
|
||||
positive_eigenvalues = np.array(eigenvalues, dtype=np.float64)
|
||||
|
||||
positive_eigenvalues[i_zero] = 1 # ones don't affect determinant
|
||||
self._log_pdet = np.sum(np.log(positive_eigenvalues), axis=-1)
|
||||
|
||||
psuedo_reciprocals = 1 / np.sqrt(positive_eigenvalues)
|
||||
psuedo_reciprocals[i_zero] = 0
|
||||
|
||||
self._LP = eigenvectors * psuedo_reciprocals
|
||||
self._LA = eigenvectors * np.sqrt(eigenvalues)
|
||||
self._rank = positive_eigenvalues.shape[-1] - i_zero.sum(axis=-1)
|
||||
self._w = eigenvalues
|
||||
self._v = eigenvectors
|
||||
self._shape = eigenvectors.shape
|
||||
self._null_basis = eigenvectors * i_zero
|
||||
# This is only used for `_support_mask`, not to decide whether
|
||||
# the covariance is singular or not.
|
||||
self._eps = _multivariate._eigvalsh_to_eps(eigenvalues) * 10**3
|
||||
self._allow_singular = True
|
||||
|
||||
def _whiten(self, x):
|
||||
return x @ self._LP
|
||||
|
||||
def _colorize(self, x):
|
||||
return x @ self._LA.T
|
||||
|
||||
@cached_property
|
||||
def _covariance(self):
|
||||
return (self._v * self._w) @ self._v.T
|
||||
|
||||
def _support_mask(self, x):
|
||||
"""
|
||||
Check whether x lies in the support of the distribution.
|
||||
"""
|
||||
residual = np.linalg.norm(x @ self._null_basis, axis=-1)
|
||||
in_support = residual < self._eps
|
||||
return in_support
|
||||
|
||||
|
||||
class CovViaPSD(Covariance):
|
||||
"""
|
||||
Representation of a covariance provided via an instance of _PSD
|
||||
"""
|
||||
|
||||
def __init__(self, psd):
|
||||
self._LP = psd.U
|
||||
self._log_pdet = psd.log_pdet
|
||||
self._rank = psd.rank
|
||||
self._covariance = psd._M
|
||||
self._shape = psd._M.shape
|
||||
self._psd = psd
|
||||
self._allow_singular = False # by default
|
||||
|
||||
def _whiten(self, x):
|
||||
return x @ self._LP
|
||||
|
||||
def _support_mask(self, x):
|
||||
return self._psd._support_mask(x)
|
||||
@@ -0,0 +1,204 @@
|
||||
import numpy as np
|
||||
from scipy.sparse import coo_matrix
|
||||
from scipy._lib._bunch import _make_tuple_bunch
|
||||
|
||||
|
||||
CrosstabResult = _make_tuple_bunch(
|
||||
"CrosstabResult", ["elements", "count"]
|
||||
)
|
||||
|
||||
|
||||
def crosstab(*args, levels=None, sparse=False):
|
||||
"""
|
||||
Return table of counts for each possible unique combination in ``*args``.
|
||||
|
||||
When ``len(args) > 1``, the array computed by this function is
|
||||
often referred to as a *contingency table* [1]_.
|
||||
|
||||
The arguments must be sequences with the same length. The second return
|
||||
value, `count`, is an integer array with ``len(args)`` dimensions. If
|
||||
`levels` is None, the shape of `count` is ``(n0, n1, ...)``, where ``nk``
|
||||
is the number of unique elements in ``args[k]``.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
*args : sequences
|
||||
A sequence of sequences whose unique aligned elements are to be
|
||||
counted. The sequences in args must all be the same length.
|
||||
levels : sequence, optional
|
||||
If `levels` is given, it must be a sequence that is the same length as
|
||||
`args`. Each element in `levels` is either a sequence or None. If it
|
||||
is a sequence, it gives the values in the corresponding sequence in
|
||||
`args` that are to be counted. If any value in the sequences in `args`
|
||||
does not occur in the corresponding sequence in `levels`, that value
|
||||
is ignored and not counted in the returned array `count`. The default
|
||||
value of `levels` for ``args[i]`` is ``np.unique(args[i])``
|
||||
sparse : bool, optional
|
||||
If True, return a sparse matrix. The matrix will be an instance of
|
||||
the `scipy.sparse.coo_matrix` class. Because SciPy's sparse matrices
|
||||
must be 2-d, only two input sequences are allowed when `sparse` is
|
||||
True. Default is False.
|
||||
|
||||
Returns
|
||||
-------
|
||||
res : CrosstabResult
|
||||
An object containing the following attributes:
|
||||
|
||||
elements : tuple of numpy.ndarrays.
|
||||
Tuple of length ``len(args)`` containing the arrays of elements
|
||||
that are counted in `count`. These can be interpreted as the
|
||||
labels of the corresponding dimensions of `count`. If `levels` was
|
||||
given, then if ``levels[i]`` is not None, ``elements[i]`` will
|
||||
hold the values given in ``levels[i]``.
|
||||
count : numpy.ndarray or scipy.sparse.coo_matrix
|
||||
Counts of the unique elements in ``zip(*args)``, stored in an
|
||||
array. Also known as a *contingency table* when ``len(args) > 1``.
|
||||
|
||||
See Also
|
||||
--------
|
||||
numpy.unique
|
||||
|
||||
Notes
|
||||
-----
|
||||
.. versionadded:: 1.7.0
|
||||
|
||||
References
|
||||
----------
|
||||
.. [1] "Contingency table", http://en.wikipedia.org/wiki/Contingency_table
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> from scipy.stats.contingency import crosstab
|
||||
|
||||
Given the lists `a` and `x`, create a contingency table that counts the
|
||||
frequencies of the corresponding pairs.
|
||||
|
||||
>>> a = ['A', 'B', 'A', 'A', 'B', 'B', 'A', 'A', 'B', 'B']
|
||||
>>> x = ['X', 'X', 'X', 'Y', 'Z', 'Z', 'Y', 'Y', 'Z', 'Z']
|
||||
>>> res = crosstab(a, x)
|
||||
>>> avals, xvals = res.elements
|
||||
>>> avals
|
||||
array(['A', 'B'], dtype='<U1')
|
||||
>>> xvals
|
||||
array(['X', 'Y', 'Z'], dtype='<U1')
|
||||
>>> res.count
|
||||
array([[2, 3, 0],
|
||||
[1, 0, 4]])
|
||||
|
||||
So `('A', 'X')` occurs twice, `('A', 'Y')` occurs three times, etc.
|
||||
|
||||
Higher dimensional contingency tables can be created.
|
||||
|
||||
>>> p = [0, 0, 0, 0, 1, 1, 1, 0, 0, 1]
|
||||
>>> res = crosstab(a, x, p)
|
||||
>>> res.count
|
||||
array([[[2, 0],
|
||||
[2, 1],
|
||||
[0, 0]],
|
||||
[[1, 0],
|
||||
[0, 0],
|
||||
[1, 3]]])
|
||||
>>> res.count.shape
|
||||
(2, 3, 2)
|
||||
|
||||
The values to be counted can be set by using the `levels` argument.
|
||||
It allows the elements of interest in each input sequence to be
|
||||
given explicitly instead finding the unique elements of the sequence.
|
||||
|
||||
For example, suppose one of the arguments is an array containing the
|
||||
answers to a survey question, with integer values 1 to 4. Even if the
|
||||
value 1 does not occur in the data, we want an entry for it in the table.
|
||||
|
||||
>>> q1 = [2, 3, 3, 2, 4, 4, 2, 3, 4, 4, 4, 3, 3, 3, 4] # 1 does not occur.
|
||||
>>> q2 = [4, 4, 2, 2, 2, 4, 1, 1, 2, 2, 4, 2, 2, 2, 4] # 3 does not occur.
|
||||
>>> options = [1, 2, 3, 4]
|
||||
>>> res = crosstab(q1, q2, levels=(options, options))
|
||||
>>> res.count
|
||||
array([[0, 0, 0, 0],
|
||||
[1, 1, 0, 1],
|
||||
[1, 4, 0, 1],
|
||||
[0, 3, 0, 3]])
|
||||
|
||||
If `levels` is given, but an element of `levels` is None, the unique values
|
||||
of the corresponding argument are used. For example,
|
||||
|
||||
>>> res = crosstab(q1, q2, levels=(None, options))
|
||||
>>> res.elements
|
||||
[array([2, 3, 4]), [1, 2, 3, 4]]
|
||||
>>> res.count
|
||||
array([[1, 1, 0, 1],
|
||||
[1, 4, 0, 1],
|
||||
[0, 3, 0, 3]])
|
||||
|
||||
If we want to ignore the pairs where 4 occurs in ``q2``, we can
|
||||
give just the values [1, 2] to `levels`, and the 4 will be ignored:
|
||||
|
||||
>>> res = crosstab(q1, q2, levels=(None, [1, 2]))
|
||||
>>> res.elements
|
||||
[array([2, 3, 4]), [1, 2]]
|
||||
>>> res.count
|
||||
array([[1, 1],
|
||||
[1, 4],
|
||||
[0, 3]])
|
||||
|
||||
Finally, let's repeat the first example, but return a sparse matrix:
|
||||
|
||||
>>> res = crosstab(a, x, sparse=True)
|
||||
>>> res.count
|
||||
<2x3 sparse matrix of type '<class 'numpy.int64'>'
|
||||
with 4 stored elements in COOrdinate format>
|
||||
>>> res.count.A
|
||||
array([[2, 3, 0],
|
||||
[1, 0, 4]])
|
||||
|
||||
"""
|
||||
nargs = len(args)
|
||||
if nargs == 0:
|
||||
raise TypeError("At least one input sequence is required.")
|
||||
|
||||
len0 = len(args[0])
|
||||
if not all(len(a) == len0 for a in args[1:]):
|
||||
raise ValueError("All input sequences must have the same length.")
|
||||
|
||||
if sparse and nargs != 2:
|
||||
raise ValueError("When `sparse` is True, only two input sequences "
|
||||
"are allowed.")
|
||||
|
||||
if levels is None:
|
||||
# Call np.unique with return_inverse=True on each argument.
|
||||
actual_levels, indices = zip(*[np.unique(a, return_inverse=True)
|
||||
for a in args])
|
||||
else:
|
||||
# `levels` is not None...
|
||||
if len(levels) != nargs:
|
||||
raise ValueError('len(levels) must equal the number of input '
|
||||
'sequences')
|
||||
|
||||
args = [np.asarray(arg) for arg in args]
|
||||
mask = np.zeros((nargs, len0), dtype=np.bool_)
|
||||
inv = np.zeros((nargs, len0), dtype=np.intp)
|
||||
actual_levels = []
|
||||
for k, (levels_list, arg) in enumerate(zip(levels, args)):
|
||||
if levels_list is None:
|
||||
levels_list, inv[k, :] = np.unique(arg, return_inverse=True)
|
||||
mask[k, :] = True
|
||||
else:
|
||||
q = arg == np.asarray(levels_list).reshape(-1, 1)
|
||||
mask[k, :] = np.any(q, axis=0)
|
||||
qnz = q.T.nonzero()
|
||||
inv[k, qnz[0]] = qnz[1]
|
||||
actual_levels.append(levels_list)
|
||||
|
||||
mask_all = mask.all(axis=0)
|
||||
indices = tuple(inv[:, mask_all])
|
||||
|
||||
if sparse:
|
||||
count = coo_matrix((np.ones(len(indices[0]), dtype=int),
|
||||
(indices[0], indices[1])))
|
||||
count.sum_duplicates()
|
||||
else:
|
||||
shape = [len(u) for u in actual_levels]
|
||||
count = np.zeros(shape, dtype=int)
|
||||
np.add.at(count, indices, 1)
|
||||
|
||||
return CrosstabResult(actual_levels, count)
|
||||
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
@@ -0,0 +1,288 @@
|
||||
"""
|
||||
Sane parameters for stats.distributions.
|
||||
"""
|
||||
import numpy as np
|
||||
|
||||
distcont = [
|
||||
['alpha', (3.5704770516650459,)],
|
||||
['anglit', ()],
|
||||
['arcsine', ()],
|
||||
['argus', (1.0,)],
|
||||
['beta', (2.3098496451481823, 0.62687954300963677)],
|
||||
['betaprime', (5, 6)],
|
||||
['bradford', (0.29891359763170633,)],
|
||||
['burr', (10.5, 4.3)],
|
||||
['burr12', (10, 4)],
|
||||
['cauchy', ()],
|
||||
['chi', (78,)],
|
||||
['chi2', (55,)],
|
||||
['cosine', ()],
|
||||
['crystalball', (2.0, 3.0)],
|
||||
['dgamma', (1.1023326088288166,)],
|
||||
['dweibull', (2.0685080649914673,)],
|
||||
['erlang', (10,)],
|
||||
['expon', ()],
|
||||
['exponnorm', (1.5,)],
|
||||
['exponpow', (2.697119160358469,)],
|
||||
['exponweib', (2.8923945291034436, 1.9505288745913174)],
|
||||
['f', (29, 18)],
|
||||
['fatiguelife', (29,)], # correction numargs = 1
|
||||
['fisk', (3.0857548622253179,)],
|
||||
['foldcauchy', (4.7164673455831894,)],
|
||||
['foldnorm', (1.9521253373555869,)],
|
||||
['gamma', (1.9932305483800778,)],
|
||||
['gausshyper', (13.763771604130699, 3.1189636648681431,
|
||||
2.5145980350183019, 5.1811649903971615)], # veryslow
|
||||
['genexpon', (9.1325976465418908, 16.231956600590632, 3.2819552690843983)],
|
||||
['genextreme', (-0.1,)],
|
||||
['gengamma', (4.4162385429431925, 3.1193091679242761)],
|
||||
['gengamma', (4.4162385429431925, -3.1193091679242761)],
|
||||
['genhalflogistic', (0.77274727809929322,)],
|
||||
['genhyperbolic', (0.5, 1.5, -0.5,)],
|
||||
['geninvgauss', (2.3, 1.5)],
|
||||
['genlogistic', (0.41192440799679475,)],
|
||||
['gennorm', (1.2988442399460265,)],
|
||||
['halfgennorm', (0.6748054997000371,)],
|
||||
['genpareto', (0.1,)], # use case with finite moments
|
||||
['gibrat', ()],
|
||||
['gompertz', (0.94743713075105251,)],
|
||||
['gumbel_l', ()],
|
||||
['gumbel_r', ()],
|
||||
['halfcauchy', ()],
|
||||
['halflogistic', ()],
|
||||
['halfnorm', ()],
|
||||
['hypsecant', ()],
|
||||
['invgamma', (4.0668996136993067,)],
|
||||
['invgauss', (0.14546264555347513,)],
|
||||
['invweibull', (10.58,)],
|
||||
['jf_skew_t', (8, 4)],
|
||||
['johnsonsb', (4.3172675099141058, 3.1837781130785063)],
|
||||
['johnsonsu', (2.554395574161155, 2.2482281679651965)],
|
||||
['kappa4', (0.0, 0.0)],
|
||||
['kappa4', (-0.1, 0.1)],
|
||||
['kappa4', (0.0, 0.1)],
|
||||
['kappa4', (0.1, 0.0)],
|
||||
['kappa3', (1.0,)],
|
||||
['ksone', (1000,)], # replace 22 by 100 to avoid failing range, ticket 956
|
||||
['kstwo', (10,)],
|
||||
['kstwobign', ()],
|
||||
['laplace', ()],
|
||||
['laplace_asymmetric', (2,)],
|
||||
['levy', ()],
|
||||
['levy_l', ()],
|
||||
['levy_stable', (1.8, -0.5)],
|
||||
['loggamma', (0.41411931826052117,)],
|
||||
['logistic', ()],
|
||||
['loglaplace', (3.2505926592051435,)],
|
||||
['lognorm', (0.95368226960575331,)],
|
||||
['loguniform', (0.01, 1.25)],
|
||||
['lomax', (1.8771398388773268,)],
|
||||
['maxwell', ()],
|
||||
['mielke', (10.4, 4.6)],
|
||||
['moyal', ()],
|
||||
['nakagami', (4.9673794866666237,)],
|
||||
['ncf', (27, 27, 0.41578441799226107)],
|
||||
['nct', (14, 0.24045031331198066)],
|
||||
['ncx2', (21, 1.0560465975116415)],
|
||||
['norm', ()],
|
||||
['norminvgauss', (1.25, 0.5)],
|
||||
['pareto', (2.621716532144454,)],
|
||||
['pearson3', (0.1,)],
|
||||
['pearson3', (-2,)],
|
||||
['powerlaw', (1.6591133289905851,)],
|
||||
['powerlaw', (0.6591133289905851,)],
|
||||
['powerlognorm', (2.1413923530064087, 0.44639540782048337)],
|
||||
['powernorm', (4.4453652254590779,)],
|
||||
['rayleigh', ()],
|
||||
['rdist', (1.6,)],
|
||||
['recipinvgauss', (0.63004267809369119,)],
|
||||
['reciprocal', (0.01, 1.25)],
|
||||
['rel_breitwigner', (36.545206797050334, )],
|
||||
['rice', (0.7749725210111873,)],
|
||||
['semicircular', ()],
|
||||
['skewcauchy', (0.5,)],
|
||||
['skewnorm', (4.0,)],
|
||||
['studentized_range', (3.0, 10.0)],
|
||||
['t', (2.7433514990818093,)],
|
||||
['trapezoid', (0.2, 0.8)],
|
||||
['triang', (0.15785029824528218,)],
|
||||
['truncexpon', (4.6907725456810478,)],
|
||||
['truncnorm', (-1.0978730080013919, 2.7306754109031979)],
|
||||
['truncnorm', (0.1, 2.)],
|
||||
['truncpareto', (1.8, 5.3)],
|
||||
['truncpareto', (2, 5)],
|
||||
['truncweibull_min', (2.5, 0.25, 1.75)],
|
||||
['tukeylambda', (3.1321477856738267,)],
|
||||
['uniform', ()],
|
||||
['vonmises', (3.9939042581071398,)],
|
||||
['vonmises_line', (3.9939042581071398,)],
|
||||
['wald', ()],
|
||||
['weibull_max', (2.8687961709100187,)],
|
||||
['weibull_min', (1.7866166930421596,)],
|
||||
['wrapcauchy', (0.031071279018614728,)]]
|
||||
|
||||
|
||||
distdiscrete = [
|
||||
['bernoulli',(0.3,)],
|
||||
['betabinom', (5, 2.3, 0.63)],
|
||||
['betanbinom', (5, 9.3, 1)],
|
||||
['binom', (5, 0.4)],
|
||||
['boltzmann',(1.4, 19)],
|
||||
['dlaplace', (0.8,)], # 0.5
|
||||
['geom', (0.5,)],
|
||||
['hypergeom',(30, 12, 6)],
|
||||
['hypergeom',(21,3,12)], # numpy.random (3,18,12) numpy ticket:921
|
||||
['hypergeom',(21,18,11)], # numpy.random (18,3,11) numpy ticket:921
|
||||
['nchypergeom_fisher', (140, 80, 60, 0.5)],
|
||||
['nchypergeom_wallenius', (140, 80, 60, 0.5)],
|
||||
['logser', (0.6,)], # re-enabled, numpy ticket:921
|
||||
['nbinom', (0.4, 0.4)], # from tickets: 583
|
||||
['nbinom', (5, 0.5)],
|
||||
['planck', (0.51,)], # 4.1
|
||||
['poisson', (0.6,)],
|
||||
['randint', (7, 31)],
|
||||
['skellam', (15, 8)],
|
||||
['zipf', (6.5,)],
|
||||
['zipfian', (0.75, 15)],
|
||||
['zipfian', (1.25, 10)],
|
||||
['yulesimon', (11.0,)],
|
||||
['nhypergeom', (20, 7, 1)]
|
||||
]
|
||||
|
||||
|
||||
invdistdiscrete = [
|
||||
# In each of the following, at least one shape parameter is invalid
|
||||
['hypergeom', (3, 3, 4)],
|
||||
['nhypergeom', (5, 2, 8)],
|
||||
['nchypergeom_fisher', (3, 3, 4, 1)],
|
||||
['nchypergeom_wallenius', (3, 3, 4, 1)],
|
||||
['bernoulli', (1.5, )],
|
||||
['binom', (10, 1.5)],
|
||||
['betabinom', (10, -0.4, -0.5)],
|
||||
['betanbinom', (10, -0.4, -0.5)],
|
||||
['boltzmann', (-1, 4)],
|
||||
['dlaplace', (-0.5, )],
|
||||
['geom', (1.5, )],
|
||||
['logser', (1.5, )],
|
||||
['nbinom', (10, 1.5)],
|
||||
['planck', (-0.5, )],
|
||||
['poisson', (-0.5, )],
|
||||
['randint', (5, 2)],
|
||||
['skellam', (-5, -2)],
|
||||
['zipf', (-2, )],
|
||||
['yulesimon', (-2, )],
|
||||
['zipfian', (-0.75, 15)]
|
||||
]
|
||||
|
||||
|
||||
invdistcont = [
|
||||
# In each of the following, at least one shape parameter is invalid
|
||||
['alpha', (-1, )],
|
||||
['anglit', ()],
|
||||
['arcsine', ()],
|
||||
['argus', (-1, )],
|
||||
['beta', (-2, 2)],
|
||||
['betaprime', (-2, 2)],
|
||||
['bradford', (-1, )],
|
||||
['burr', (-1, 1)],
|
||||
['burr12', (-1, 1)],
|
||||
['cauchy', ()],
|
||||
['chi', (-1, )],
|
||||
['chi2', (-1, )],
|
||||
['cosine', ()],
|
||||
['crystalball', (-1, 2)],
|
||||
['dgamma', (-1, )],
|
||||
['dweibull', (-1, )],
|
||||
['erlang', (-1, )],
|
||||
['expon', ()],
|
||||
['exponnorm', (-1, )],
|
||||
['exponweib', (1, -1)],
|
||||
['exponpow', (-1, )],
|
||||
['f', (10, -10)],
|
||||
['fatiguelife', (-1, )],
|
||||
['fisk', (-1, )],
|
||||
['foldcauchy', (-1, )],
|
||||
['foldnorm', (-1, )],
|
||||
['genlogistic', (-1, )],
|
||||
['gennorm', (-1, )],
|
||||
['genpareto', (np.inf, )],
|
||||
['genexpon', (1, 2, -3)],
|
||||
['genextreme', (np.inf, )],
|
||||
['genhyperbolic', (0.5, -0.5, -1.5,)],
|
||||
['gausshyper', (1, 2, 3, -4)],
|
||||
['gamma', (-1, )],
|
||||
['gengamma', (-1, 0)],
|
||||
['genhalflogistic', (-1, )],
|
||||
['geninvgauss', (1, 0)],
|
||||
['gibrat', ()],
|
||||
['gompertz', (-1, )],
|
||||
['gumbel_r', ()],
|
||||
['gumbel_l', ()],
|
||||
['halfcauchy', ()],
|
||||
['halflogistic', ()],
|
||||
['halfnorm', ()],
|
||||
['halfgennorm', (-1, )],
|
||||
['hypsecant', ()],
|
||||
['invgamma', (-1, )],
|
||||
['invgauss', (-1, )],
|
||||
['invweibull', (-1, )],
|
||||
['jf_skew_t', (-1, 0)],
|
||||
['johnsonsb', (1, -2)],
|
||||
['johnsonsu', (1, -2)],
|
||||
['kappa4', (np.nan, 0)],
|
||||
['kappa3', (-1, )],
|
||||
['ksone', (-1, )],
|
||||
['kstwo', (-1, )],
|
||||
['kstwobign', ()],
|
||||
['laplace', ()],
|
||||
['laplace_asymmetric', (-1, )],
|
||||
['levy', ()],
|
||||
['levy_l', ()],
|
||||
['levy_stable', (-1, 1)],
|
||||
['logistic', ()],
|
||||
['loggamma', (-1, )],
|
||||
['loglaplace', (-1, )],
|
||||
['lognorm', (-1, )],
|
||||
['loguniform', (10, 5)],
|
||||
['lomax', (-1, )],
|
||||
['maxwell', ()],
|
||||
['mielke', (1, -2)],
|
||||
['moyal', ()],
|
||||
['nakagami', (-1, )],
|
||||
['ncx2', (-1, 2)],
|
||||
['ncf', (10, 20, -1)],
|
||||
['nct', (-1, 2)],
|
||||
['norm', ()],
|
||||
['norminvgauss', (5, -10)],
|
||||
['pareto', (-1, )],
|
||||
['pearson3', (np.nan, )],
|
||||
['powerlaw', (-1, )],
|
||||
['powerlognorm', (1, -2)],
|
||||
['powernorm', (-1, )],
|
||||
['rdist', (-1, )],
|
||||
['rayleigh', ()],
|
||||
['rice', (-1, )],
|
||||
['recipinvgauss', (-1, )],
|
||||
['semicircular', ()],
|
||||
['skewnorm', (np.inf, )],
|
||||
['studentized_range', (-1, 1)],
|
||||
['rel_breitwigner', (-2, )],
|
||||
['t', (-1, )],
|
||||
['trapezoid', (0, 2)],
|
||||
['triang', (2, )],
|
||||
['truncexpon', (-1, )],
|
||||
['truncnorm', (10, 5)],
|
||||
['truncpareto', (-1, 5)],
|
||||
['truncpareto', (1.8, .5)],
|
||||
['truncweibull_min', (-2.5, 0.25, 1.75)],
|
||||
['tukeylambda', (np.nan, )],
|
||||
['uniform', ()],
|
||||
['vonmises', (-1, )],
|
||||
['vonmises_line', (-1, )],
|
||||
['wald', ()],
|
||||
['weibull_min', (-1, )],
|
||||
['weibull_max', (-1, )],
|
||||
['wrapcauchy', (2, )],
|
||||
['reciprocal', (15, 10)],
|
||||
['skewcauchy', (2, )]
|
||||
]
|
||||
@@ -0,0 +1,423 @@
|
||||
"""
|
||||
Created on Fri Apr 2 09:06:05 2021
|
||||
|
||||
@author: matth
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
import math
|
||||
import numpy as np
|
||||
from scipy import special
|
||||
from ._axis_nan_policy import _axis_nan_policy_factory, _broadcast_arrays
|
||||
|
||||
__all__ = ['entropy', 'differential_entropy']
|
||||
|
||||
|
||||
@_axis_nan_policy_factory(
|
||||
lambda x: x,
|
||||
n_samples=lambda kwgs: (
|
||||
2 if ("qk" in kwgs and kwgs["qk"] is not None)
|
||||
else 1
|
||||
),
|
||||
n_outputs=1, result_to_tuple=lambda x: (x,), paired=True,
|
||||
too_small=-1 # entropy doesn't have too small inputs
|
||||
)
|
||||
def entropy(pk: np.typing.ArrayLike,
|
||||
qk: np.typing.ArrayLike | None = None,
|
||||
base: float | None = None,
|
||||
axis: int = 0
|
||||
) -> np.number | np.ndarray:
|
||||
"""
|
||||
Calculate the Shannon entropy/relative entropy of given distribution(s).
|
||||
|
||||
If only probabilities `pk` are given, the Shannon entropy is calculated as
|
||||
``H = -sum(pk * log(pk))``.
|
||||
|
||||
If `qk` is not None, then compute the relative entropy
|
||||
``D = sum(pk * log(pk / qk))``. This quantity is also known
|
||||
as the Kullback-Leibler divergence.
|
||||
|
||||
This routine will normalize `pk` and `qk` if they don't sum to 1.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
pk : array_like
|
||||
Defines the (discrete) distribution. Along each axis-slice of ``pk``,
|
||||
element ``i`` is the (possibly unnormalized) probability of event
|
||||
``i``.
|
||||
qk : array_like, optional
|
||||
Sequence against which the relative entropy is computed. Should be in
|
||||
the same format as `pk`.
|
||||
base : float, optional
|
||||
The logarithmic base to use, defaults to ``e`` (natural logarithm).
|
||||
axis : int, optional
|
||||
The axis along which the entropy is calculated. Default is 0.
|
||||
|
||||
Returns
|
||||
-------
|
||||
S : {float, array_like}
|
||||
The calculated entropy.
|
||||
|
||||
Notes
|
||||
-----
|
||||
Informally, the Shannon entropy quantifies the expected uncertainty
|
||||
inherent in the possible outcomes of a discrete random variable.
|
||||
For example,
|
||||
if messages consisting of sequences of symbols from a set are to be
|
||||
encoded and transmitted over a noiseless channel, then the Shannon entropy
|
||||
``H(pk)`` gives a tight lower bound for the average number of units of
|
||||
information needed per symbol if the symbols occur with frequencies
|
||||
governed by the discrete distribution `pk` [1]_. The choice of base
|
||||
determines the choice of units; e.g., ``e`` for nats, ``2`` for bits, etc.
|
||||
|
||||
The relative entropy, ``D(pk|qk)``, quantifies the increase in the average
|
||||
number of units of information needed per symbol if the encoding is
|
||||
optimized for the probability distribution `qk` instead of the true
|
||||
distribution `pk`. Informally, the relative entropy quantifies the expected
|
||||
excess in surprise experienced if one believes the true distribution is
|
||||
`qk` when it is actually `pk`.
|
||||
|
||||
A related quantity, the cross entropy ``CE(pk, qk)``, satisfies the
|
||||
equation ``CE(pk, qk) = H(pk) + D(pk|qk)`` and can also be calculated with
|
||||
the formula ``CE = -sum(pk * log(qk))``. It gives the average
|
||||
number of units of information needed per symbol if an encoding is
|
||||
optimized for the probability distribution `qk` when the true distribution
|
||||
is `pk`. It is not computed directly by `entropy`, but it can be computed
|
||||
using two calls to the function (see Examples).
|
||||
|
||||
See [2]_ for more information.
|
||||
|
||||
References
|
||||
----------
|
||||
.. [1] Shannon, C.E. (1948), A Mathematical Theory of Communication.
|
||||
Bell System Technical Journal, 27: 379-423.
|
||||
https://doi.org/10.1002/j.1538-7305.1948.tb01338.x
|
||||
.. [2] Thomas M. Cover and Joy A. Thomas. 2006. Elements of Information
|
||||
Theory (Wiley Series in Telecommunications and Signal Processing).
|
||||
Wiley-Interscience, USA.
|
||||
|
||||
|
||||
Examples
|
||||
--------
|
||||
The outcome of a fair coin is the most uncertain:
|
||||
|
||||
>>> import numpy as np
|
||||
>>> from scipy.stats import entropy
|
||||
>>> base = 2 # work in units of bits
|
||||
>>> pk = np.array([1/2, 1/2]) # fair coin
|
||||
>>> H = entropy(pk, base=base)
|
||||
>>> H
|
||||
1.0
|
||||
>>> H == -np.sum(pk * np.log(pk)) / np.log(base)
|
||||
True
|
||||
|
||||
The outcome of a biased coin is less uncertain:
|
||||
|
||||
>>> qk = np.array([9/10, 1/10]) # biased coin
|
||||
>>> entropy(qk, base=base)
|
||||
0.46899559358928117
|
||||
|
||||
The relative entropy between the fair coin and biased coin is calculated
|
||||
as:
|
||||
|
||||
>>> D = entropy(pk, qk, base=base)
|
||||
>>> D
|
||||
0.7369655941662062
|
||||
>>> D == np.sum(pk * np.log(pk/qk)) / np.log(base)
|
||||
True
|
||||
|
||||
The cross entropy can be calculated as the sum of the entropy and
|
||||
relative entropy`:
|
||||
|
||||
>>> CE = entropy(pk, base=base) + entropy(pk, qk, base=base)
|
||||
>>> CE
|
||||
1.736965594166206
|
||||
>>> CE == -np.sum(pk * np.log(qk)) / np.log(base)
|
||||
True
|
||||
|
||||
"""
|
||||
if base is not None and base <= 0:
|
||||
raise ValueError("`base` must be a positive number or `None`.")
|
||||
|
||||
pk = np.asarray(pk)
|
||||
with np.errstate(invalid='ignore'):
|
||||
pk = 1.0*pk / np.sum(pk, axis=axis, keepdims=True)
|
||||
if qk is None:
|
||||
vec = special.entr(pk)
|
||||
else:
|
||||
qk = np.asarray(qk)
|
||||
pk, qk = _broadcast_arrays((pk, qk), axis=None) # don't ignore any axes
|
||||
sum_kwargs = dict(axis=axis, keepdims=True)
|
||||
qk = 1.0*qk / np.sum(qk, **sum_kwargs) # type: ignore[operator, call-overload]
|
||||
vec = special.rel_entr(pk, qk)
|
||||
S = np.sum(vec, axis=axis)
|
||||
if base is not None:
|
||||
S /= np.log(base)
|
||||
return S
|
||||
|
||||
|
||||
def _differential_entropy_is_too_small(samples, kwargs, axis=-1):
|
||||
values = samples[0]
|
||||
n = values.shape[axis]
|
||||
window_length = kwargs.get("window_length",
|
||||
math.floor(math.sqrt(n) + 0.5))
|
||||
if not 2 <= 2 * window_length < n:
|
||||
return True
|
||||
return False
|
||||
|
||||
|
||||
@_axis_nan_policy_factory(
|
||||
lambda x: x, n_outputs=1, result_to_tuple=lambda x: (x,),
|
||||
too_small=_differential_entropy_is_too_small
|
||||
)
|
||||
def differential_entropy(
|
||||
values: np.typing.ArrayLike,
|
||||
*,
|
||||
window_length: int | None = None,
|
||||
base: float | None = None,
|
||||
axis: int = 0,
|
||||
method: str = "auto",
|
||||
) -> np.number | np.ndarray:
|
||||
r"""Given a sample of a distribution, estimate the differential entropy.
|
||||
|
||||
Several estimation methods are available using the `method` parameter. By
|
||||
default, a method is selected based the size of the sample.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
values : sequence
|
||||
Sample from a continuous distribution.
|
||||
window_length : int, optional
|
||||
Window length for computing Vasicek estimate. Must be an integer
|
||||
between 1 and half of the sample size. If ``None`` (the default), it
|
||||
uses the heuristic value
|
||||
|
||||
.. math::
|
||||
\left \lfloor \sqrt{n} + 0.5 \right \rfloor
|
||||
|
||||
where :math:`n` is the sample size. This heuristic was originally
|
||||
proposed in [2]_ and has become common in the literature.
|
||||
base : float, optional
|
||||
The logarithmic base to use, defaults to ``e`` (natural logarithm).
|
||||
axis : int, optional
|
||||
The axis along which the differential entropy is calculated.
|
||||
Default is 0.
|
||||
method : {'vasicek', 'van es', 'ebrahimi', 'correa', 'auto'}, optional
|
||||
The method used to estimate the differential entropy from the sample.
|
||||
Default is ``'auto'``. See Notes for more information.
|
||||
|
||||
Returns
|
||||
-------
|
||||
entropy : float
|
||||
The calculated differential entropy.
|
||||
|
||||
Notes
|
||||
-----
|
||||
This function will converge to the true differential entropy in the limit
|
||||
|
||||
.. math::
|
||||
n \to \infty, \quad m \to \infty, \quad \frac{m}{n} \to 0
|
||||
|
||||
The optimal choice of ``window_length`` for a given sample size depends on
|
||||
the (unknown) distribution. Typically, the smoother the density of the
|
||||
distribution, the larger the optimal value of ``window_length`` [1]_.
|
||||
|
||||
The following options are available for the `method` parameter.
|
||||
|
||||
* ``'vasicek'`` uses the estimator presented in [1]_. This is
|
||||
one of the first and most influential estimators of differential entropy.
|
||||
* ``'van es'`` uses the bias-corrected estimator presented in [3]_, which
|
||||
is not only consistent but, under some conditions, asymptotically normal.
|
||||
* ``'ebrahimi'`` uses an estimator presented in [4]_, which was shown
|
||||
in simulation to have smaller bias and mean squared error than
|
||||
the Vasicek estimator.
|
||||
* ``'correa'`` uses the estimator presented in [5]_ based on local linear
|
||||
regression. In a simulation study, it had consistently smaller mean
|
||||
square error than the Vasiceck estimator, but it is more expensive to
|
||||
compute.
|
||||
* ``'auto'`` selects the method automatically (default). Currently,
|
||||
this selects ``'van es'`` for very small samples (<10), ``'ebrahimi'``
|
||||
for moderate sample sizes (11-1000), and ``'vasicek'`` for larger
|
||||
samples, but this behavior is subject to change in future versions.
|
||||
|
||||
All estimators are implemented as described in [6]_.
|
||||
|
||||
References
|
||||
----------
|
||||
.. [1] Vasicek, O. (1976). A test for normality based on sample entropy.
|
||||
Journal of the Royal Statistical Society:
|
||||
Series B (Methodological), 38(1), 54-59.
|
||||
.. [2] Crzcgorzewski, P., & Wirczorkowski, R. (1999). Entropy-based
|
||||
goodness-of-fit test for exponentiality. Communications in
|
||||
Statistics-Theory and Methods, 28(5), 1183-1202.
|
||||
.. [3] Van Es, B. (1992). Estimating functionals related to a density by a
|
||||
class of statistics based on spacings. Scandinavian Journal of
|
||||
Statistics, 61-72.
|
||||
.. [4] Ebrahimi, N., Pflughoeft, K., & Soofi, E. S. (1994). Two measures
|
||||
of sample entropy. Statistics & Probability Letters, 20(3), 225-234.
|
||||
.. [5] Correa, J. C. (1995). A new estimator of entropy. Communications
|
||||
in Statistics-Theory and Methods, 24(10), 2439-2449.
|
||||
.. [6] Noughabi, H. A. (2015). Entropy Estimation Using Numerical Methods.
|
||||
Annals of Data Science, 2(2), 231-241.
|
||||
https://link.springer.com/article/10.1007/s40745-015-0045-9
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> import numpy as np
|
||||
>>> from scipy.stats import differential_entropy, norm
|
||||
|
||||
Entropy of a standard normal distribution:
|
||||
|
||||
>>> rng = np.random.default_rng()
|
||||
>>> values = rng.standard_normal(100)
|
||||
>>> differential_entropy(values)
|
||||
1.3407817436640392
|
||||
|
||||
Compare with the true entropy:
|
||||
|
||||
>>> float(norm.entropy())
|
||||
1.4189385332046727
|
||||
|
||||
For several sample sizes between 5 and 1000, compare the accuracy of
|
||||
the ``'vasicek'``, ``'van es'``, and ``'ebrahimi'`` methods. Specifically,
|
||||
compare the root mean squared error (over 1000 trials) between the estimate
|
||||
and the true differential entropy of the distribution.
|
||||
|
||||
>>> from scipy import stats
|
||||
>>> import matplotlib.pyplot as plt
|
||||
>>>
|
||||
>>>
|
||||
>>> def rmse(res, expected):
|
||||
... '''Root mean squared error'''
|
||||
... return np.sqrt(np.mean((res - expected)**2))
|
||||
>>>
|
||||
>>>
|
||||
>>> a, b = np.log10(5), np.log10(1000)
|
||||
>>> ns = np.round(np.logspace(a, b, 10)).astype(int)
|
||||
>>> reps = 1000 # number of repetitions for each sample size
|
||||
>>> expected = stats.expon.entropy()
|
||||
>>>
|
||||
>>> method_errors = {'vasicek': [], 'van es': [], 'ebrahimi': []}
|
||||
>>> for method in method_errors:
|
||||
... for n in ns:
|
||||
... rvs = stats.expon.rvs(size=(reps, n), random_state=rng)
|
||||
... res = stats.differential_entropy(rvs, method=method, axis=-1)
|
||||
... error = rmse(res, expected)
|
||||
... method_errors[method].append(error)
|
||||
>>>
|
||||
>>> for method, errors in method_errors.items():
|
||||
... plt.loglog(ns, errors, label=method)
|
||||
>>>
|
||||
>>> plt.legend()
|
||||
>>> plt.xlabel('sample size')
|
||||
>>> plt.ylabel('RMSE (1000 trials)')
|
||||
>>> plt.title('Entropy Estimator Error (Exponential Distribution)')
|
||||
|
||||
"""
|
||||
values = np.asarray(values)
|
||||
values = np.moveaxis(values, axis, -1)
|
||||
n = values.shape[-1] # number of observations
|
||||
|
||||
if window_length is None:
|
||||
window_length = math.floor(math.sqrt(n) + 0.5)
|
||||
|
||||
if not 2 <= 2 * window_length < n:
|
||||
raise ValueError(
|
||||
f"Window length ({window_length}) must be positive and less "
|
||||
f"than half the sample size ({n}).",
|
||||
)
|
||||
|
||||
if base is not None and base <= 0:
|
||||
raise ValueError("`base` must be a positive number or `None`.")
|
||||
|
||||
sorted_data = np.sort(values, axis=-1)
|
||||
|
||||
methods = {"vasicek": _vasicek_entropy,
|
||||
"van es": _van_es_entropy,
|
||||
"correa": _correa_entropy,
|
||||
"ebrahimi": _ebrahimi_entropy,
|
||||
"auto": _vasicek_entropy}
|
||||
method = method.lower()
|
||||
if method not in methods:
|
||||
message = f"`method` must be one of {set(methods)}"
|
||||
raise ValueError(message)
|
||||
|
||||
if method == "auto":
|
||||
if n <= 10:
|
||||
method = 'van es'
|
||||
elif n <= 1000:
|
||||
method = 'ebrahimi'
|
||||
else:
|
||||
method = 'vasicek'
|
||||
|
||||
res = methods[method](sorted_data, window_length)
|
||||
|
||||
if base is not None:
|
||||
res /= np.log(base)
|
||||
|
||||
return res
|
||||
|
||||
|
||||
def _pad_along_last_axis(X, m):
|
||||
"""Pad the data for computing the rolling window difference."""
|
||||
# scales a bit better than method in _vasicek_like_entropy
|
||||
shape = np.array(X.shape)
|
||||
shape[-1] = m
|
||||
Xl = np.broadcast_to(X[..., [0]], shape) # [0] vs 0 to maintain shape
|
||||
Xr = np.broadcast_to(X[..., [-1]], shape)
|
||||
return np.concatenate((Xl, X, Xr), axis=-1)
|
||||
|
||||
|
||||
def _vasicek_entropy(X, m):
|
||||
"""Compute the Vasicek estimator as described in [6] Eq. 1.3."""
|
||||
n = X.shape[-1]
|
||||
X = _pad_along_last_axis(X, m)
|
||||
differences = X[..., 2 * m:] - X[..., : -2 * m:]
|
||||
logs = np.log(n/(2*m) * differences)
|
||||
return np.mean(logs, axis=-1)
|
||||
|
||||
|
||||
def _van_es_entropy(X, m):
|
||||
"""Compute the van Es estimator as described in [6]."""
|
||||
# No equation number, but referred to as HVE_mn.
|
||||
# Typo: there should be a log within the summation.
|
||||
n = X.shape[-1]
|
||||
difference = X[..., m:] - X[..., :-m]
|
||||
term1 = 1/(n-m) * np.sum(np.log((n+1)/m * difference), axis=-1)
|
||||
k = np.arange(m, n+1)
|
||||
return term1 + np.sum(1/k) + np.log(m) - np.log(n+1)
|
||||
|
||||
|
||||
def _ebrahimi_entropy(X, m):
|
||||
"""Compute the Ebrahimi estimator as described in [6]."""
|
||||
# No equation number, but referred to as HE_mn
|
||||
n = X.shape[-1]
|
||||
X = _pad_along_last_axis(X, m)
|
||||
|
||||
differences = X[..., 2 * m:] - X[..., : -2 * m:]
|
||||
|
||||
i = np.arange(1, n+1).astype(float)
|
||||
ci = np.ones_like(i)*2
|
||||
ci[i <= m] = 1 + (i[i <= m] - 1)/m
|
||||
ci[i >= n - m + 1] = 1 + (n - i[i >= n-m+1])/m
|
||||
|
||||
logs = np.log(n * differences / (ci * m))
|
||||
return np.mean(logs, axis=-1)
|
||||
|
||||
|
||||
def _correa_entropy(X, m):
|
||||
"""Compute the Correa estimator as described in [6]."""
|
||||
# No equation number, but referred to as HC_mn
|
||||
n = X.shape[-1]
|
||||
X = _pad_along_last_axis(X, m)
|
||||
|
||||
i = np.arange(1, n+1)
|
||||
dj = np.arange(-m, m+1)[:, None]
|
||||
j = i + dj
|
||||
j0 = j + m - 1 # 0-indexed version of j
|
||||
|
||||
Xibar = np.mean(X[..., j0], axis=-2, keepdims=True)
|
||||
difference = X[..., j0] - Xibar
|
||||
num = np.sum(difference*dj, axis=-2) # dj is d-i
|
||||
den = n*np.sum(difference**2, axis=-2)
|
||||
return -np.mean(np.log(num/den), axis=-1)
|
||||
File diff suppressed because it is too large
Load Diff
@@ -0,0 +1,27 @@
|
||||
import pathlib
|
||||
import subprocess
|
||||
import sys
|
||||
import os
|
||||
import argparse
|
||||
|
||||
|
||||
def make_boost(outdir):
|
||||
# Call code generator inside _boost directory
|
||||
code_gen = pathlib.Path(__file__).parent / '_boost/include/code_gen.py'
|
||||
subprocess.run([sys.executable, str(code_gen), '-o', outdir],
|
||||
check=True)
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
parser = argparse.ArgumentParser()
|
||||
parser.add_argument("-o", "--outdir", type=str,
|
||||
help="Path to the output directory")
|
||||
args = parser.parse_args()
|
||||
|
||||
if not args.outdir:
|
||||
raise ValueError("A path to the output directory is required")
|
||||
else:
|
||||
# Meson build
|
||||
srcdir_abs = pathlib.Path(os.path.abspath(os.path.dirname(__file__)))
|
||||
outdir_abs = pathlib.Path(os.getcwd()) / args.outdir
|
||||
make_boost(outdir_abs)
|
||||
File diff suppressed because it is too large
Load Diff
@@ -0,0 +1,728 @@
|
||||
#-------------------------------------------------------------------------------
|
||||
#
|
||||
# Define classes for (uni/multi)-variate kernel density estimation.
|
||||
#
|
||||
# Currently, only Gaussian kernels are implemented.
|
||||
#
|
||||
# Written by: Robert Kern
|
||||
#
|
||||
# Date: 2004-08-09
|
||||
#
|
||||
# Modified: 2005-02-10 by Robert Kern.
|
||||
# Contributed to SciPy
|
||||
# 2005-10-07 by Robert Kern.
|
||||
# Some fixes to match the new scipy_core
|
||||
#
|
||||
# Copyright 2004-2005 by Enthought, Inc.
|
||||
#
|
||||
#-------------------------------------------------------------------------------
|
||||
|
||||
# Standard library imports.
|
||||
import warnings
|
||||
|
||||
# SciPy imports.
|
||||
from scipy import linalg, special
|
||||
from scipy._lib._util import check_random_state
|
||||
|
||||
from numpy import (asarray, atleast_2d, reshape, zeros, newaxis, exp, pi,
|
||||
sqrt, ravel, power, atleast_1d, squeeze, sum, transpose,
|
||||
ones, cov)
|
||||
import numpy as np
|
||||
|
||||
# Local imports.
|
||||
from . import _mvn
|
||||
from ._stats import gaussian_kernel_estimate, gaussian_kernel_estimate_log
|
||||
|
||||
# deprecated import to be removed in SciPy 1.13.0
|
||||
from scipy.special import logsumexp # noqa: F401
|
||||
|
||||
|
||||
__all__ = ['gaussian_kde']
|
||||
|
||||
|
||||
class gaussian_kde:
|
||||
"""Representation of a kernel-density estimate using Gaussian kernels.
|
||||
|
||||
Kernel density estimation is a way to estimate the probability density
|
||||
function (PDF) of a random variable in a non-parametric way.
|
||||
`gaussian_kde` works for both uni-variate and multi-variate data. It
|
||||
includes automatic bandwidth determination. The estimation works best for
|
||||
a unimodal distribution; bimodal or multi-modal distributions tend to be
|
||||
oversmoothed.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
dataset : array_like
|
||||
Datapoints to estimate from. In case of univariate data this is a 1-D
|
||||
array, otherwise a 2-D array with shape (# of dims, # of data).
|
||||
bw_method : str, scalar or callable, optional
|
||||
The method used to calculate the estimator bandwidth. This can be
|
||||
'scott', 'silverman', a scalar constant or a callable. If a scalar,
|
||||
this will be used directly as `kde.factor`. If a callable, it should
|
||||
take a `gaussian_kde` instance as only parameter and return a scalar.
|
||||
If None (default), 'scott' is used. See Notes for more details.
|
||||
weights : array_like, optional
|
||||
weights of datapoints. This must be the same shape as dataset.
|
||||
If None (default), the samples are assumed to be equally weighted
|
||||
|
||||
Attributes
|
||||
----------
|
||||
dataset : ndarray
|
||||
The dataset with which `gaussian_kde` was initialized.
|
||||
d : int
|
||||
Number of dimensions.
|
||||
n : int
|
||||
Number of datapoints.
|
||||
neff : int
|
||||
Effective number of datapoints.
|
||||
|
||||
.. versionadded:: 1.2.0
|
||||
factor : float
|
||||
The bandwidth factor, obtained from `kde.covariance_factor`. The square
|
||||
of `kde.factor` multiplies the covariance matrix of the data in the kde
|
||||
estimation.
|
||||
covariance : ndarray
|
||||
The covariance matrix of `dataset`, scaled by the calculated bandwidth
|
||||
(`kde.factor`).
|
||||
inv_cov : ndarray
|
||||
The inverse of `covariance`.
|
||||
|
||||
Methods
|
||||
-------
|
||||
evaluate
|
||||
__call__
|
||||
integrate_gaussian
|
||||
integrate_box_1d
|
||||
integrate_box
|
||||
integrate_kde
|
||||
pdf
|
||||
logpdf
|
||||
resample
|
||||
set_bandwidth
|
||||
covariance_factor
|
||||
|
||||
Notes
|
||||
-----
|
||||
Bandwidth selection strongly influences the estimate obtained from the KDE
|
||||
(much more so than the actual shape of the kernel). Bandwidth selection
|
||||
can be done by a "rule of thumb", by cross-validation, by "plug-in
|
||||
methods" or by other means; see [3]_, [4]_ for reviews. `gaussian_kde`
|
||||
uses a rule of thumb, the default is Scott's Rule.
|
||||
|
||||
Scott's Rule [1]_, implemented as `scotts_factor`, is::
|
||||
|
||||
n**(-1./(d+4)),
|
||||
|
||||
with ``n`` the number of data points and ``d`` the number of dimensions.
|
||||
In the case of unequally weighted points, `scotts_factor` becomes::
|
||||
|
||||
neff**(-1./(d+4)),
|
||||
|
||||
with ``neff`` the effective number of datapoints.
|
||||
Silverman's Rule [2]_, implemented as `silverman_factor`, is::
|
||||
|
||||
(n * (d + 2) / 4.)**(-1. / (d + 4)).
|
||||
|
||||
or in the case of unequally weighted points::
|
||||
|
||||
(neff * (d + 2) / 4.)**(-1. / (d + 4)).
|
||||
|
||||
Good general descriptions of kernel density estimation can be found in [1]_
|
||||
and [2]_, the mathematics for this multi-dimensional implementation can be
|
||||
found in [1]_.
|
||||
|
||||
With a set of weighted samples, the effective number of datapoints ``neff``
|
||||
is defined by::
|
||||
|
||||
neff = sum(weights)^2 / sum(weights^2)
|
||||
|
||||
as detailed in [5]_.
|
||||
|
||||
`gaussian_kde` does not currently support data that lies in a
|
||||
lower-dimensional subspace of the space in which it is expressed. For such
|
||||
data, consider performing principle component analysis / dimensionality
|
||||
reduction and using `gaussian_kde` with the transformed data.
|
||||
|
||||
References
|
||||
----------
|
||||
.. [1] D.W. Scott, "Multivariate Density Estimation: Theory, Practice, and
|
||||
Visualization", John Wiley & Sons, New York, Chicester, 1992.
|
||||
.. [2] B.W. Silverman, "Density Estimation for Statistics and Data
|
||||
Analysis", Vol. 26, Monographs on Statistics and Applied Probability,
|
||||
Chapman and Hall, London, 1986.
|
||||
.. [3] B.A. Turlach, "Bandwidth Selection in Kernel Density Estimation: A
|
||||
Review", CORE and Institut de Statistique, Vol. 19, pp. 1-33, 1993.
|
||||
.. [4] D.M. Bashtannyk and R.J. Hyndman, "Bandwidth selection for kernel
|
||||
conditional density estimation", Computational Statistics & Data
|
||||
Analysis, Vol. 36, pp. 279-298, 2001.
|
||||
.. [5] Gray P. G., 1969, Journal of the Royal Statistical Society.
|
||||
Series A (General), 132, 272
|
||||
|
||||
Examples
|
||||
--------
|
||||
Generate some random two-dimensional data:
|
||||
|
||||
>>> import numpy as np
|
||||
>>> from scipy import stats
|
||||
>>> def measure(n):
|
||||
... "Measurement model, return two coupled measurements."
|
||||
... m1 = np.random.normal(size=n)
|
||||
... m2 = np.random.normal(scale=0.5, size=n)
|
||||
... return m1+m2, m1-m2
|
||||
|
||||
>>> m1, m2 = measure(2000)
|
||||
>>> xmin = m1.min()
|
||||
>>> xmax = m1.max()
|
||||
>>> ymin = m2.min()
|
||||
>>> ymax = m2.max()
|
||||
|
||||
Perform a kernel density estimate on the data:
|
||||
|
||||
>>> X, Y = np.mgrid[xmin:xmax:100j, ymin:ymax:100j]
|
||||
>>> positions = np.vstack([X.ravel(), Y.ravel()])
|
||||
>>> values = np.vstack([m1, m2])
|
||||
>>> kernel = stats.gaussian_kde(values)
|
||||
>>> Z = np.reshape(kernel(positions).T, X.shape)
|
||||
|
||||
Plot the results:
|
||||
|
||||
>>> import matplotlib.pyplot as plt
|
||||
>>> fig, ax = plt.subplots()
|
||||
>>> ax.imshow(np.rot90(Z), cmap=plt.cm.gist_earth_r,
|
||||
... extent=[xmin, xmax, ymin, ymax])
|
||||
>>> ax.plot(m1, m2, 'k.', markersize=2)
|
||||
>>> ax.set_xlim([xmin, xmax])
|
||||
>>> ax.set_ylim([ymin, ymax])
|
||||
>>> plt.show()
|
||||
|
||||
"""
|
||||
def __init__(self, dataset, bw_method=None, weights=None):
|
||||
self.dataset = atleast_2d(asarray(dataset))
|
||||
if not self.dataset.size > 1:
|
||||
raise ValueError("`dataset` input should have multiple elements.")
|
||||
|
||||
self.d, self.n = self.dataset.shape
|
||||
|
||||
if weights is not None:
|
||||
self._weights = atleast_1d(weights).astype(float)
|
||||
self._weights /= sum(self._weights)
|
||||
if self.weights.ndim != 1:
|
||||
raise ValueError("`weights` input should be one-dimensional.")
|
||||
if len(self._weights) != self.n:
|
||||
raise ValueError("`weights` input should be of length n")
|
||||
self._neff = 1/sum(self._weights**2)
|
||||
|
||||
# This can be converted to a warning once gh-10205 is resolved
|
||||
if self.d > self.n:
|
||||
msg = ("Number of dimensions is greater than number of samples. "
|
||||
"This results in a singular data covariance matrix, which "
|
||||
"cannot be treated using the algorithms implemented in "
|
||||
"`gaussian_kde`. Note that `gaussian_kde` interprets each "
|
||||
"*column* of `dataset` to be a point; consider transposing "
|
||||
"the input to `dataset`.")
|
||||
raise ValueError(msg)
|
||||
|
||||
try:
|
||||
self.set_bandwidth(bw_method=bw_method)
|
||||
except linalg.LinAlgError as e:
|
||||
msg = ("The data appears to lie in a lower-dimensional subspace "
|
||||
"of the space in which it is expressed. This has resulted "
|
||||
"in a singular data covariance matrix, which cannot be "
|
||||
"treated using the algorithms implemented in "
|
||||
"`gaussian_kde`. Consider performing principle component "
|
||||
"analysis / dimensionality reduction and using "
|
||||
"`gaussian_kde` with the transformed data.")
|
||||
raise linalg.LinAlgError(msg) from e
|
||||
|
||||
def evaluate(self, points):
|
||||
"""Evaluate the estimated pdf on a set of points.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
points : (# of dimensions, # of points)-array
|
||||
Alternatively, a (# of dimensions,) vector can be passed in and
|
||||
treated as a single point.
|
||||
|
||||
Returns
|
||||
-------
|
||||
values : (# of points,)-array
|
||||
The values at each point.
|
||||
|
||||
Raises
|
||||
------
|
||||
ValueError : if the dimensionality of the input points is different than
|
||||
the dimensionality of the KDE.
|
||||
|
||||
"""
|
||||
points = atleast_2d(asarray(points))
|
||||
|
||||
d, m = points.shape
|
||||
if d != self.d:
|
||||
if d == 1 and m == self.d:
|
||||
# points was passed in as a row vector
|
||||
points = reshape(points, (self.d, 1))
|
||||
m = 1
|
||||
else:
|
||||
msg = (f"points have dimension {d}, "
|
||||
f"dataset has dimension {self.d}")
|
||||
raise ValueError(msg)
|
||||
|
||||
output_dtype, spec = _get_output_dtype(self.covariance, points)
|
||||
result = gaussian_kernel_estimate[spec](
|
||||
self.dataset.T, self.weights[:, None],
|
||||
points.T, self.cho_cov, output_dtype)
|
||||
|
||||
return result[:, 0]
|
||||
|
||||
__call__ = evaluate
|
||||
|
||||
def integrate_gaussian(self, mean, cov):
|
||||
"""
|
||||
Multiply estimated density by a multivariate Gaussian and integrate
|
||||
over the whole space.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
mean : aray_like
|
||||
A 1-D array, specifying the mean of the Gaussian.
|
||||
cov : array_like
|
||||
A 2-D array, specifying the covariance matrix of the Gaussian.
|
||||
|
||||
Returns
|
||||
-------
|
||||
result : scalar
|
||||
The value of the integral.
|
||||
|
||||
Raises
|
||||
------
|
||||
ValueError
|
||||
If the mean or covariance of the input Gaussian differs from
|
||||
the KDE's dimensionality.
|
||||
|
||||
"""
|
||||
mean = atleast_1d(squeeze(mean))
|
||||
cov = atleast_2d(cov)
|
||||
|
||||
if mean.shape != (self.d,):
|
||||
raise ValueError("mean does not have dimension %s" % self.d)
|
||||
if cov.shape != (self.d, self.d):
|
||||
raise ValueError("covariance does not have dimension %s" % self.d)
|
||||
|
||||
# make mean a column vector
|
||||
mean = mean[:, newaxis]
|
||||
|
||||
sum_cov = self.covariance + cov
|
||||
|
||||
# This will raise LinAlgError if the new cov matrix is not s.p.d
|
||||
# cho_factor returns (ndarray, bool) where bool is a flag for whether
|
||||
# or not ndarray is upper or lower triangular
|
||||
sum_cov_chol = linalg.cho_factor(sum_cov)
|
||||
|
||||
diff = self.dataset - mean
|
||||
tdiff = linalg.cho_solve(sum_cov_chol, diff)
|
||||
|
||||
sqrt_det = np.prod(np.diagonal(sum_cov_chol[0]))
|
||||
norm_const = power(2 * pi, sum_cov.shape[0] / 2.0) * sqrt_det
|
||||
|
||||
energies = sum(diff * tdiff, axis=0) / 2.0
|
||||
result = sum(exp(-energies)*self.weights, axis=0) / norm_const
|
||||
|
||||
return result
|
||||
|
||||
def integrate_box_1d(self, low, high):
|
||||
"""
|
||||
Computes the integral of a 1D pdf between two bounds.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
low : scalar
|
||||
Lower bound of integration.
|
||||
high : scalar
|
||||
Upper bound of integration.
|
||||
|
||||
Returns
|
||||
-------
|
||||
value : scalar
|
||||
The result of the integral.
|
||||
|
||||
Raises
|
||||
------
|
||||
ValueError
|
||||
If the KDE is over more than one dimension.
|
||||
|
||||
"""
|
||||
if self.d != 1:
|
||||
raise ValueError("integrate_box_1d() only handles 1D pdfs")
|
||||
|
||||
stdev = ravel(sqrt(self.covariance))[0]
|
||||
|
||||
normalized_low = ravel((low - self.dataset) / stdev)
|
||||
normalized_high = ravel((high - self.dataset) / stdev)
|
||||
|
||||
value = np.sum(self.weights*(
|
||||
special.ndtr(normalized_high) -
|
||||
special.ndtr(normalized_low)))
|
||||
return value
|
||||
|
||||
def integrate_box(self, low_bounds, high_bounds, maxpts=None):
|
||||
"""Computes the integral of a pdf over a rectangular interval.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
low_bounds : array_like
|
||||
A 1-D array containing the lower bounds of integration.
|
||||
high_bounds : array_like
|
||||
A 1-D array containing the upper bounds of integration.
|
||||
maxpts : int, optional
|
||||
The maximum number of points to use for integration.
|
||||
|
||||
Returns
|
||||
-------
|
||||
value : scalar
|
||||
The result of the integral.
|
||||
|
||||
"""
|
||||
if maxpts is not None:
|
||||
extra_kwds = {'maxpts': maxpts}
|
||||
else:
|
||||
extra_kwds = {}
|
||||
|
||||
value, inform = _mvn.mvnun_weighted(low_bounds, high_bounds,
|
||||
self.dataset, self.weights,
|
||||
self.covariance, **extra_kwds)
|
||||
if inform:
|
||||
msg = ('An integral in _mvn.mvnun requires more points than %s' %
|
||||
(self.d * 1000))
|
||||
warnings.warn(msg, stacklevel=2)
|
||||
|
||||
return value
|
||||
|
||||
def integrate_kde(self, other):
|
||||
"""
|
||||
Computes the integral of the product of this kernel density estimate
|
||||
with another.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
other : gaussian_kde instance
|
||||
The other kde.
|
||||
|
||||
Returns
|
||||
-------
|
||||
value : scalar
|
||||
The result of the integral.
|
||||
|
||||
Raises
|
||||
------
|
||||
ValueError
|
||||
If the KDEs have different dimensionality.
|
||||
|
||||
"""
|
||||
if other.d != self.d:
|
||||
raise ValueError("KDEs are not the same dimensionality")
|
||||
|
||||
# we want to iterate over the smallest number of points
|
||||
if other.n < self.n:
|
||||
small = other
|
||||
large = self
|
||||
else:
|
||||
small = self
|
||||
large = other
|
||||
|
||||
sum_cov = small.covariance + large.covariance
|
||||
sum_cov_chol = linalg.cho_factor(sum_cov)
|
||||
result = 0.0
|
||||
for i in range(small.n):
|
||||
mean = small.dataset[:, i, newaxis]
|
||||
diff = large.dataset - mean
|
||||
tdiff = linalg.cho_solve(sum_cov_chol, diff)
|
||||
|
||||
energies = sum(diff * tdiff, axis=0) / 2.0
|
||||
result += sum(exp(-energies)*large.weights, axis=0)*small.weights[i]
|
||||
|
||||
sqrt_det = np.prod(np.diagonal(sum_cov_chol[0]))
|
||||
norm_const = power(2 * pi, sum_cov.shape[0] / 2.0) * sqrt_det
|
||||
|
||||
result /= norm_const
|
||||
|
||||
return result
|
||||
|
||||
def resample(self, size=None, seed=None):
|
||||
"""Randomly sample a dataset from the estimated pdf.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
size : int, optional
|
||||
The number of samples to draw. If not provided, then the size is
|
||||
the same as the effective number of samples in the underlying
|
||||
dataset.
|
||||
seed : {None, int, `numpy.random.Generator`, `numpy.random.RandomState`}, optional
|
||||
If `seed` is None (or `np.random`), the `numpy.random.RandomState`
|
||||
singleton is used.
|
||||
If `seed` is an int, a new ``RandomState`` instance is used,
|
||||
seeded with `seed`.
|
||||
If `seed` is already a ``Generator`` or ``RandomState`` instance then
|
||||
that instance is used.
|
||||
|
||||
Returns
|
||||
-------
|
||||
resample : (self.d, `size`) ndarray
|
||||
The sampled dataset.
|
||||
|
||||
""" # numpy/numpydoc#87 # noqa: E501
|
||||
if size is None:
|
||||
size = int(self.neff)
|
||||
|
||||
random_state = check_random_state(seed)
|
||||
norm = transpose(random_state.multivariate_normal(
|
||||
zeros((self.d,), float), self.covariance, size=size
|
||||
))
|
||||
indices = random_state.choice(self.n, size=size, p=self.weights)
|
||||
means = self.dataset[:, indices]
|
||||
|
||||
return means + norm
|
||||
|
||||
def scotts_factor(self):
|
||||
"""Compute Scott's factor.
|
||||
|
||||
Returns
|
||||
-------
|
||||
s : float
|
||||
Scott's factor.
|
||||
"""
|
||||
return power(self.neff, -1./(self.d+4))
|
||||
|
||||
def silverman_factor(self):
|
||||
"""Compute the Silverman factor.
|
||||
|
||||
Returns
|
||||
-------
|
||||
s : float
|
||||
The silverman factor.
|
||||
"""
|
||||
return power(self.neff*(self.d+2.0)/4.0, -1./(self.d+4))
|
||||
|
||||
# Default method to calculate bandwidth, can be overwritten by subclass
|
||||
covariance_factor = scotts_factor
|
||||
covariance_factor.__doc__ = """Computes the coefficient (`kde.factor`) that
|
||||
multiplies the data covariance matrix to obtain the kernel covariance
|
||||
matrix. The default is `scotts_factor`. A subclass can overwrite this
|
||||
method to provide a different method, or set it through a call to
|
||||
`kde.set_bandwidth`."""
|
||||
|
||||
def set_bandwidth(self, bw_method=None):
|
||||
"""Compute the estimator bandwidth with given method.
|
||||
|
||||
The new bandwidth calculated after a call to `set_bandwidth` is used
|
||||
for subsequent evaluations of the estimated density.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
bw_method : str, scalar or callable, optional
|
||||
The method used to calculate the estimator bandwidth. This can be
|
||||
'scott', 'silverman', a scalar constant or a callable. If a
|
||||
scalar, this will be used directly as `kde.factor`. If a callable,
|
||||
it should take a `gaussian_kde` instance as only parameter and
|
||||
return a scalar. If None (default), nothing happens; the current
|
||||
`kde.covariance_factor` method is kept.
|
||||
|
||||
Notes
|
||||
-----
|
||||
.. versionadded:: 0.11
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> import numpy as np
|
||||
>>> import scipy.stats as stats
|
||||
>>> x1 = np.array([-7, -5, 1, 4, 5.])
|
||||
>>> kde = stats.gaussian_kde(x1)
|
||||
>>> xs = np.linspace(-10, 10, num=50)
|
||||
>>> y1 = kde(xs)
|
||||
>>> kde.set_bandwidth(bw_method='silverman')
|
||||
>>> y2 = kde(xs)
|
||||
>>> kde.set_bandwidth(bw_method=kde.factor / 3.)
|
||||
>>> y3 = kde(xs)
|
||||
|
||||
>>> import matplotlib.pyplot as plt
|
||||
>>> fig, ax = plt.subplots()
|
||||
>>> ax.plot(x1, np.full(x1.shape, 1 / (4. * x1.size)), 'bo',
|
||||
... label='Data points (rescaled)')
|
||||
>>> ax.plot(xs, y1, label='Scott (default)')
|
||||
>>> ax.plot(xs, y2, label='Silverman')
|
||||
>>> ax.plot(xs, y3, label='Const (1/3 * Silverman)')
|
||||
>>> ax.legend()
|
||||
>>> plt.show()
|
||||
|
||||
"""
|
||||
if bw_method is None:
|
||||
pass
|
||||
elif bw_method == 'scott':
|
||||
self.covariance_factor = self.scotts_factor
|
||||
elif bw_method == 'silverman':
|
||||
self.covariance_factor = self.silverman_factor
|
||||
elif np.isscalar(bw_method) and not isinstance(bw_method, str):
|
||||
self._bw_method = 'use constant'
|
||||
self.covariance_factor = lambda: bw_method
|
||||
elif callable(bw_method):
|
||||
self._bw_method = bw_method
|
||||
self.covariance_factor = lambda: self._bw_method(self)
|
||||
else:
|
||||
msg = "`bw_method` should be 'scott', 'silverman', a scalar " \
|
||||
"or a callable."
|
||||
raise ValueError(msg)
|
||||
|
||||
self._compute_covariance()
|
||||
|
||||
def _compute_covariance(self):
|
||||
"""Computes the covariance matrix for each Gaussian kernel using
|
||||
covariance_factor().
|
||||
"""
|
||||
self.factor = self.covariance_factor()
|
||||
# Cache covariance and Cholesky decomp of covariance
|
||||
if not hasattr(self, '_data_cho_cov'):
|
||||
self._data_covariance = atleast_2d(cov(self.dataset, rowvar=1,
|
||||
bias=False,
|
||||
aweights=self.weights))
|
||||
self._data_cho_cov = linalg.cholesky(self._data_covariance,
|
||||
lower=True)
|
||||
|
||||
self.covariance = self._data_covariance * self.factor**2
|
||||
self.cho_cov = (self._data_cho_cov * self.factor).astype(np.float64)
|
||||
self.log_det = 2*np.log(np.diag(self.cho_cov
|
||||
* np.sqrt(2*pi))).sum()
|
||||
|
||||
@property
|
||||
def inv_cov(self):
|
||||
# Re-compute from scratch each time because I'm not sure how this is
|
||||
# used in the wild. (Perhaps users change the `dataset`, since it's
|
||||
# not a private attribute?) `_compute_covariance` used to recalculate
|
||||
# all these, so we'll recalculate everything now that this is a
|
||||
# a property.
|
||||
self.factor = self.covariance_factor()
|
||||
self._data_covariance = atleast_2d(cov(self.dataset, rowvar=1,
|
||||
bias=False, aweights=self.weights))
|
||||
return linalg.inv(self._data_covariance) / self.factor**2
|
||||
|
||||
def pdf(self, x):
|
||||
"""
|
||||
Evaluate the estimated pdf on a provided set of points.
|
||||
|
||||
Notes
|
||||
-----
|
||||
This is an alias for `gaussian_kde.evaluate`. See the ``evaluate``
|
||||
docstring for more details.
|
||||
|
||||
"""
|
||||
return self.evaluate(x)
|
||||
|
||||
def logpdf(self, x):
|
||||
"""
|
||||
Evaluate the log of the estimated pdf on a provided set of points.
|
||||
"""
|
||||
points = atleast_2d(x)
|
||||
|
||||
d, m = points.shape
|
||||
if d != self.d:
|
||||
if d == 1 and m == self.d:
|
||||
# points was passed in as a row vector
|
||||
points = reshape(points, (self.d, 1))
|
||||
m = 1
|
||||
else:
|
||||
msg = (f"points have dimension {d}, "
|
||||
f"dataset has dimension {self.d}")
|
||||
raise ValueError(msg)
|
||||
|
||||
output_dtype, spec = _get_output_dtype(self.covariance, points)
|
||||
result = gaussian_kernel_estimate_log[spec](
|
||||
self.dataset.T, self.weights[:, None],
|
||||
points.T, self.cho_cov, output_dtype)
|
||||
|
||||
return result[:, 0]
|
||||
|
||||
def marginal(self, dimensions):
|
||||
"""Return a marginal KDE distribution
|
||||
|
||||
Parameters
|
||||
----------
|
||||
dimensions : int or 1-d array_like
|
||||
The dimensions of the multivariate distribution corresponding
|
||||
with the marginal variables, that is, the indices of the dimensions
|
||||
that are being retained. The other dimensions are marginalized out.
|
||||
|
||||
Returns
|
||||
-------
|
||||
marginal_kde : gaussian_kde
|
||||
An object representing the marginal distribution.
|
||||
|
||||
Notes
|
||||
-----
|
||||
.. versionadded:: 1.10.0
|
||||
|
||||
"""
|
||||
|
||||
dims = np.atleast_1d(dimensions)
|
||||
|
||||
if not np.issubdtype(dims.dtype, np.integer):
|
||||
msg = ("Elements of `dimensions` must be integers - the indices "
|
||||
"of the marginal variables being retained.")
|
||||
raise ValueError(msg)
|
||||
|
||||
n = len(self.dataset) # number of dimensions
|
||||
original_dims = dims.copy()
|
||||
|
||||
dims[dims < 0] = n + dims[dims < 0]
|
||||
|
||||
if len(np.unique(dims)) != len(dims):
|
||||
msg = ("All elements of `dimensions` must be unique.")
|
||||
raise ValueError(msg)
|
||||
|
||||
i_invalid = (dims < 0) | (dims >= n)
|
||||
if np.any(i_invalid):
|
||||
msg = (f"Dimensions {original_dims[i_invalid]} are invalid "
|
||||
f"for a distribution in {n} dimensions.")
|
||||
raise ValueError(msg)
|
||||
|
||||
dataset = self.dataset[dims]
|
||||
weights = self.weights
|
||||
|
||||
return gaussian_kde(dataset, bw_method=self.covariance_factor(),
|
||||
weights=weights)
|
||||
|
||||
@property
|
||||
def weights(self):
|
||||
try:
|
||||
return self._weights
|
||||
except AttributeError:
|
||||
self._weights = ones(self.n)/self.n
|
||||
return self._weights
|
||||
|
||||
@property
|
||||
def neff(self):
|
||||
try:
|
||||
return self._neff
|
||||
except AttributeError:
|
||||
self._neff = 1/sum(self.weights**2)
|
||||
return self._neff
|
||||
|
||||
|
||||
def _get_output_dtype(covariance, points):
|
||||
"""
|
||||
Calculates the output dtype and the "spec" (=C type name).
|
||||
|
||||
This was necessary in order to deal with the fused types in the Cython
|
||||
routine `gaussian_kernel_estimate`. See gh-10824 for details.
|
||||
"""
|
||||
output_dtype = np.common_type(covariance, points)
|
||||
itemsize = np.dtype(output_dtype).itemsize
|
||||
if itemsize == 4:
|
||||
spec = 'float'
|
||||
elif itemsize == 8:
|
||||
spec = 'double'
|
||||
elif itemsize in (12, 16):
|
||||
spec = 'long double'
|
||||
else:
|
||||
raise ValueError(
|
||||
f"{output_dtype} has unexpected item size: {itemsize}"
|
||||
)
|
||||
|
||||
return output_dtype, spec
|
||||
@@ -0,0 +1,600 @@
|
||||
# Compute the two-sided one-sample Kolmogorov-Smirnov Prob(Dn <= d) where:
|
||||
# D_n = sup_x{|F_n(x) - F(x)|},
|
||||
# F_n(x) is the empirical CDF for a sample of size n {x_i: i=1,...,n},
|
||||
# F(x) is the CDF of a probability distribution.
|
||||
#
|
||||
# Exact methods:
|
||||
# Prob(D_n >= d) can be computed via a matrix algorithm of Durbin[1]
|
||||
# or a recursion algorithm due to Pomeranz[2].
|
||||
# Marsaglia, Tsang & Wang[3] gave a computation-efficient way to perform
|
||||
# the Durbin algorithm.
|
||||
# D_n >= d <==> D_n+ >= d or D_n- >= d (the one-sided K-S statistics), hence
|
||||
# Prob(D_n >= d) = 2*Prob(D_n+ >= d) - Prob(D_n+ >= d and D_n- >= d).
|
||||
# For d > 0.5, the latter intersection probability is 0.
|
||||
#
|
||||
# Approximate methods:
|
||||
# For d close to 0.5, ignoring that intersection term may still give a
|
||||
# reasonable approximation.
|
||||
# Li-Chien[4] and Korolyuk[5] gave an asymptotic formula extending
|
||||
# Kolmogorov's initial asymptotic, suitable for large d. (See
|
||||
# scipy.special.kolmogorov for that asymptotic)
|
||||
# Pelz-Good[6] used the functional equation for Jacobi theta functions to
|
||||
# transform the Li-Chien/Korolyuk formula produce a computational formula
|
||||
# suitable for small d.
|
||||
#
|
||||
# Simard and L'Ecuyer[7] provided an algorithm to decide when to use each of
|
||||
# the above approaches and it is that which is used here.
|
||||
#
|
||||
# Other approaches:
|
||||
# Carvalho[8] optimizes Durbin's matrix algorithm for large values of d.
|
||||
# Moscovich and Nadler[9] use FFTs to compute the convolutions.
|
||||
|
||||
# References:
|
||||
# [1] Durbin J (1968).
|
||||
# "The Probability that the Sample Distribution Function Lies Between Two
|
||||
# Parallel Straight Lines."
|
||||
# Annals of Mathematical Statistics, 39, 398-411.
|
||||
# [2] Pomeranz J (1974).
|
||||
# "Exact Cumulative Distribution of the Kolmogorov-Smirnov Statistic for
|
||||
# Small Samples (Algorithm 487)."
|
||||
# Communications of the ACM, 17(12), 703-704.
|
||||
# [3] Marsaglia G, Tsang WW, Wang J (2003).
|
||||
# "Evaluating Kolmogorov's Distribution."
|
||||
# Journal of Statistical Software, 8(18), 1-4.
|
||||
# [4] LI-CHIEN, C. (1956).
|
||||
# "On the exact distribution of the statistics of A. N. Kolmogorov and
|
||||
# their asymptotic expansion."
|
||||
# Acta Matematica Sinica, 6, 55-81.
|
||||
# [5] KOROLYUK, V. S. (1960).
|
||||
# "Asymptotic analysis of the distribution of the maximum deviation in
|
||||
# the Bernoulli scheme."
|
||||
# Theor. Probability Appl., 4, 339-366.
|
||||
# [6] Pelz W, Good IJ (1976).
|
||||
# "Approximating the Lower Tail-areas of the Kolmogorov-Smirnov One-sample
|
||||
# Statistic."
|
||||
# Journal of the Royal Statistical Society, Series B, 38(2), 152-156.
|
||||
# [7] Simard, R., L'Ecuyer, P. (2011)
|
||||
# "Computing the Two-Sided Kolmogorov-Smirnov Distribution",
|
||||
# Journal of Statistical Software, Vol 39, 11, 1-18.
|
||||
# [8] Carvalho, Luis (2015)
|
||||
# "An Improved Evaluation of Kolmogorov's Distribution"
|
||||
# Journal of Statistical Software, Code Snippets; Vol 65(3), 1-8.
|
||||
# [9] Amit Moscovich, Boaz Nadler (2017)
|
||||
# "Fast calculation of boundary crossing probabilities for Poisson
|
||||
# processes",
|
||||
# Statistics & Probability Letters, Vol 123, 177-182.
|
||||
|
||||
|
||||
import numpy as np
|
||||
import scipy.special
|
||||
import scipy.special._ufuncs as scu
|
||||
from scipy._lib._finite_differences import _derivative
|
||||
|
||||
_E128 = 128
|
||||
_EP128 = np.ldexp(np.longdouble(1), _E128)
|
||||
_EM128 = np.ldexp(np.longdouble(1), -_E128)
|
||||
|
||||
_SQRT2PI = np.sqrt(2 * np.pi)
|
||||
_LOG_2PI = np.log(2 * np.pi)
|
||||
_MIN_LOG = -708
|
||||
_SQRT3 = np.sqrt(3)
|
||||
_PI_SQUARED = np.pi ** 2
|
||||
_PI_FOUR = np.pi ** 4
|
||||
_PI_SIX = np.pi ** 6
|
||||
|
||||
# [Lifted from _loggamma.pxd.] If B_m are the Bernoulli numbers,
|
||||
# then Stirling coeffs are B_{2j}/(2j)/(2j-1) for j=8,...1.
|
||||
_STIRLING_COEFFS = [-2.955065359477124183e-2, 6.4102564102564102564e-3,
|
||||
-1.9175269175269175269e-3, 8.4175084175084175084e-4,
|
||||
-5.952380952380952381e-4, 7.9365079365079365079e-4,
|
||||
-2.7777777777777777778e-3, 8.3333333333333333333e-2]
|
||||
|
||||
|
||||
def _log_nfactorial_div_n_pow_n(n):
|
||||
# Computes n! / n**n
|
||||
# = (n-1)! / n**(n-1)
|
||||
# Uses Stirling's approximation, but removes n*log(n) up-front to
|
||||
# avoid subtractive cancellation.
|
||||
# = log(n)/2 - n + log(sqrt(2pi)) + sum B_{2j}/(2j)/(2j-1)/n**(2j-1)
|
||||
rn = 1.0/n
|
||||
return np.log(n)/2 - n + _LOG_2PI/2 + rn * np.polyval(_STIRLING_COEFFS, rn/n)
|
||||
|
||||
|
||||
def _clip_prob(p):
|
||||
"""clips a probability to range 0<=p<=1."""
|
||||
return np.clip(p, 0.0, 1.0)
|
||||
|
||||
|
||||
def _select_and_clip_prob(cdfprob, sfprob, cdf=True):
|
||||
"""Selects either the CDF or SF, and then clips to range 0<=p<=1."""
|
||||
p = np.where(cdf, cdfprob, sfprob)
|
||||
return _clip_prob(p)
|
||||
|
||||
|
||||
def _kolmogn_DMTW(n, d, cdf=True):
|
||||
r"""Computes the Kolmogorov CDF: Pr(D_n <= d) using the MTW approach to
|
||||
the Durbin matrix algorithm.
|
||||
|
||||
Durbin (1968); Marsaglia, Tsang, Wang (2003). [1], [3].
|
||||
"""
|
||||
# Write d = (k-h)/n, where k is positive integer and 0 <= h < 1
|
||||
# Generate initial matrix H of size m*m where m=(2k-1)
|
||||
# Compute k-th row of (n!/n^n) * H^n, scaling intermediate results.
|
||||
# Requires memory O(m^2) and computation O(m^2 log(n)).
|
||||
# Most suitable for small m.
|
||||
|
||||
if d >= 1.0:
|
||||
return _select_and_clip_prob(1.0, 0.0, cdf)
|
||||
nd = n * d
|
||||
if nd <= 0.5:
|
||||
return _select_and_clip_prob(0.0, 1.0, cdf)
|
||||
k = int(np.ceil(nd))
|
||||
h = k - nd
|
||||
m = 2 * k - 1
|
||||
|
||||
H = np.zeros([m, m])
|
||||
|
||||
# Initialize: v is first column (and last row) of H
|
||||
# v[j] = (1-h^(j+1)/(j+1)! (except for v[-1])
|
||||
# w[j] = 1/(j)!
|
||||
# q = k-th row of H (actually i!/n^i*H^i)
|
||||
intm = np.arange(1, m + 1)
|
||||
v = 1.0 - h ** intm
|
||||
w = np.empty(m)
|
||||
fac = 1.0
|
||||
for j in intm:
|
||||
w[j - 1] = fac
|
||||
fac /= j # This might underflow. Isn't a problem.
|
||||
v[j - 1] *= fac
|
||||
tt = max(2 * h - 1.0, 0)**m - 2*h**m
|
||||
v[-1] = (1.0 + tt) * fac
|
||||
|
||||
for i in range(1, m):
|
||||
H[i - 1:, i] = w[:m - i + 1]
|
||||
H[:, 0] = v
|
||||
H[-1, :] = np.flip(v, axis=0)
|
||||
|
||||
Hpwr = np.eye(np.shape(H)[0]) # Holds intermediate powers of H
|
||||
nn = n
|
||||
expnt = 0 # Scaling of Hpwr
|
||||
Hexpnt = 0 # Scaling of H
|
||||
while nn > 0:
|
||||
if nn % 2:
|
||||
Hpwr = np.matmul(Hpwr, H)
|
||||
expnt += Hexpnt
|
||||
H = np.matmul(H, H)
|
||||
Hexpnt *= 2
|
||||
# Scale as needed.
|
||||
if np.abs(H[k - 1, k - 1]) > _EP128:
|
||||
H /= _EP128
|
||||
Hexpnt += _E128
|
||||
nn = nn // 2
|
||||
|
||||
p = Hpwr[k - 1, k - 1]
|
||||
|
||||
# Multiply by n!/n^n
|
||||
for i in range(1, n + 1):
|
||||
p = i * p / n
|
||||
if np.abs(p) < _EM128:
|
||||
p *= _EP128
|
||||
expnt -= _E128
|
||||
|
||||
# unscale
|
||||
if expnt != 0:
|
||||
p = np.ldexp(p, expnt)
|
||||
|
||||
return _select_and_clip_prob(p, 1.0-p, cdf)
|
||||
|
||||
|
||||
def _pomeranz_compute_j1j2(i, n, ll, ceilf, roundf):
|
||||
"""Compute the endpoints of the interval for row i."""
|
||||
if i == 0:
|
||||
j1, j2 = -ll - ceilf - 1, ll + ceilf - 1
|
||||
else:
|
||||
# i + 1 = 2*ip1div2 + ip1mod2
|
||||
ip1div2, ip1mod2 = divmod(i + 1, 2)
|
||||
if ip1mod2 == 0: # i is odd
|
||||
if ip1div2 == n + 1:
|
||||
j1, j2 = n - ll - ceilf - 1, n + ll + ceilf - 1
|
||||
else:
|
||||
j1, j2 = ip1div2 - 1 - ll - roundf - 1, ip1div2 + ll - 1 + ceilf - 1
|
||||
else:
|
||||
j1, j2 = ip1div2 - 1 - ll - 1, ip1div2 + ll + roundf - 1
|
||||
|
||||
return max(j1 + 2, 0), min(j2, n)
|
||||
|
||||
|
||||
def _kolmogn_Pomeranz(n, x, cdf=True):
|
||||
r"""Computes Pr(D_n <= d) using the Pomeranz recursion algorithm.
|
||||
|
||||
Pomeranz (1974) [2]
|
||||
"""
|
||||
|
||||
# V is n*(2n+2) matrix.
|
||||
# Each row is convolution of the previous row and probabilities from a
|
||||
# Poisson distribution.
|
||||
# Desired CDF probability is n! V[n-1, 2n+1] (final entry in final row).
|
||||
# Only two rows are needed at any given stage:
|
||||
# - Call them V0 and V1.
|
||||
# - Swap each iteration
|
||||
# Only a few (contiguous) entries in each row can be non-zero.
|
||||
# - Keep track of start and end (j1 and j2 below)
|
||||
# - V0s and V1s track the start in the two rows
|
||||
# Scale intermediate results as needed.
|
||||
# Only a few different Poisson distributions can occur
|
||||
t = n * x
|
||||
ll = int(np.floor(t))
|
||||
f = 1.0 * (t - ll) # fractional part of t
|
||||
g = min(f, 1.0 - f)
|
||||
ceilf = (1 if f > 0 else 0)
|
||||
roundf = (1 if f > 0.5 else 0)
|
||||
npwrs = 2 * (ll + 1) # Maximum number of powers needed in convolutions
|
||||
gpower = np.empty(npwrs) # gpower = (g/n)^m/m!
|
||||
twogpower = np.empty(npwrs) # twogpower = (2g/n)^m/m!
|
||||
onem2gpower = np.empty(npwrs) # onem2gpower = ((1-2g)/n)^m/m!
|
||||
# gpower etc are *almost* Poisson probs, just missing normalizing factor.
|
||||
|
||||
gpower[0] = 1.0
|
||||
twogpower[0] = 1.0
|
||||
onem2gpower[0] = 1.0
|
||||
expnt = 0
|
||||
g_over_n, two_g_over_n, one_minus_two_g_over_n = g/n, 2*g/n, (1 - 2*g)/n
|
||||
for m in range(1, npwrs):
|
||||
gpower[m] = gpower[m - 1] * g_over_n / m
|
||||
twogpower[m] = twogpower[m - 1] * two_g_over_n / m
|
||||
onem2gpower[m] = onem2gpower[m - 1] * one_minus_two_g_over_n / m
|
||||
|
||||
V0 = np.zeros([npwrs])
|
||||
V1 = np.zeros([npwrs])
|
||||
V1[0] = 1 # first row
|
||||
V0s, V1s = 0, 0 # start indices of the two rows
|
||||
|
||||
j1, j2 = _pomeranz_compute_j1j2(0, n, ll, ceilf, roundf)
|
||||
for i in range(1, 2 * n + 2):
|
||||
# Preserve j1, V1, V1s, V0s from last iteration
|
||||
k1 = j1
|
||||
V0, V1 = V1, V0
|
||||
V0s, V1s = V1s, V0s
|
||||
V1.fill(0.0)
|
||||
j1, j2 = _pomeranz_compute_j1j2(i, n, ll, ceilf, roundf)
|
||||
if i == 1 or i == 2 * n + 1:
|
||||
pwrs = gpower
|
||||
else:
|
||||
pwrs = (twogpower if i % 2 else onem2gpower)
|
||||
ln2 = j2 - k1 + 1
|
||||
if ln2 > 0:
|
||||
conv = np.convolve(V0[k1 - V0s:k1 - V0s + ln2], pwrs[:ln2])
|
||||
conv_start = j1 - k1 # First index to use from conv
|
||||
conv_len = j2 - j1 + 1 # Number of entries to use from conv
|
||||
V1[:conv_len] = conv[conv_start:conv_start + conv_len]
|
||||
# Scale to avoid underflow.
|
||||
if 0 < np.max(V1) < _EM128:
|
||||
V1 *= _EP128
|
||||
expnt -= _E128
|
||||
V1s = V0s + j1 - k1
|
||||
|
||||
# multiply by n!
|
||||
ans = V1[n - V1s]
|
||||
for m in range(1, n + 1):
|
||||
if np.abs(ans) > _EP128:
|
||||
ans *= _EM128
|
||||
expnt += _E128
|
||||
ans *= m
|
||||
|
||||
# Undo any intermediate scaling
|
||||
if expnt != 0:
|
||||
ans = np.ldexp(ans, expnt)
|
||||
ans = _select_and_clip_prob(ans, 1.0 - ans, cdf)
|
||||
return ans
|
||||
|
||||
|
||||
def _kolmogn_PelzGood(n, x, cdf=True):
|
||||
"""Computes the Pelz-Good approximation to Prob(Dn <= x) with 0<=x<=1.
|
||||
|
||||
Start with Li-Chien, Korolyuk approximation:
|
||||
Prob(Dn <= x) ~ K0(z) + K1(z)/sqrt(n) + K2(z)/n + K3(z)/n**1.5
|
||||
where z = x*sqrt(n).
|
||||
Transform each K_(z) using Jacobi theta functions into a form suitable
|
||||
for small z.
|
||||
Pelz-Good (1976). [6]
|
||||
"""
|
||||
if x <= 0.0:
|
||||
return _select_and_clip_prob(0.0, 1.0, cdf=cdf)
|
||||
if x >= 1.0:
|
||||
return _select_and_clip_prob(1.0, 0.0, cdf=cdf)
|
||||
|
||||
z = np.sqrt(n) * x
|
||||
zsquared, zthree, zfour, zsix = z**2, z**3, z**4, z**6
|
||||
|
||||
qlog = -_PI_SQUARED / 8 / zsquared
|
||||
if qlog < _MIN_LOG: # z ~ 0.041743441416853426
|
||||
return _select_and_clip_prob(0.0, 1.0, cdf=cdf)
|
||||
|
||||
q = np.exp(qlog)
|
||||
|
||||
# Coefficients of terms in the sums for K1, K2 and K3
|
||||
k1a = -zsquared
|
||||
k1b = _PI_SQUARED / 4
|
||||
|
||||
k2a = 6 * zsix + 2 * zfour
|
||||
k2b = (2 * zfour - 5 * zsquared) * _PI_SQUARED / 4
|
||||
k2c = _PI_FOUR * (1 - 2 * zsquared) / 16
|
||||
|
||||
k3d = _PI_SIX * (5 - 30 * zsquared) / 64
|
||||
k3c = _PI_FOUR * (-60 * zsquared + 212 * zfour) / 16
|
||||
k3b = _PI_SQUARED * (135 * zfour - 96 * zsix) / 4
|
||||
k3a = -30 * zsix - 90 * z**8
|
||||
|
||||
K0to3 = np.zeros(4)
|
||||
# Use a Horner scheme to evaluate sum c_i q^(i^2)
|
||||
# Reduces to a sum over odd integers.
|
||||
maxk = int(np.ceil(16 * z / np.pi))
|
||||
for k in range(maxk, 0, -1):
|
||||
m = 2 * k - 1
|
||||
msquared, mfour, msix = m**2, m**4, m**6
|
||||
qpower = np.power(q, 8 * k)
|
||||
coeffs = np.array([1.0,
|
||||
k1a + k1b*msquared,
|
||||
k2a + k2b*msquared + k2c*mfour,
|
||||
k3a + k3b*msquared + k3c*mfour + k3d*msix])
|
||||
K0to3 *= qpower
|
||||
K0to3 += coeffs
|
||||
K0to3 *= q
|
||||
K0to3 *= _SQRT2PI
|
||||
# z**10 > 0 as z > 0.04
|
||||
K0to3 /= np.array([z, 6 * zfour, 72 * z**7, 6480 * z**10])
|
||||
|
||||
# Now do the other sum over the other terms, all integers k
|
||||
# K_2: (pi^2 k^2) q^(k^2),
|
||||
# K_3: (3pi^2 k^2 z^2 - pi^4 k^4)*q^(k^2)
|
||||
# Don't expect much subtractive cancellation so use direct calculation
|
||||
q = np.exp(-_PI_SQUARED / 2 / zsquared)
|
||||
ks = np.arange(maxk, 0, -1)
|
||||
ksquared = ks ** 2
|
||||
sqrt3z = _SQRT3 * z
|
||||
kspi = np.pi * ks
|
||||
qpwers = q ** ksquared
|
||||
k2extra = np.sum(ksquared * qpwers)
|
||||
k2extra *= _PI_SQUARED * _SQRT2PI/(-36 * zthree)
|
||||
K0to3[2] += k2extra
|
||||
k3extra = np.sum((sqrt3z + kspi) * (sqrt3z - kspi) * ksquared * qpwers)
|
||||
k3extra *= _PI_SQUARED * _SQRT2PI/(216 * zsix)
|
||||
K0to3[3] += k3extra
|
||||
powers_of_n = np.power(n * 1.0, np.arange(len(K0to3)) / 2.0)
|
||||
K0to3 /= powers_of_n
|
||||
|
||||
if not cdf:
|
||||
K0to3 *= -1
|
||||
K0to3[0] += 1
|
||||
|
||||
Ksum = sum(K0to3)
|
||||
return Ksum
|
||||
|
||||
|
||||
def _kolmogn(n, x, cdf=True):
|
||||
"""Computes the CDF(or SF) for the two-sided Kolmogorov-Smirnov statistic.
|
||||
|
||||
x must be of type float, n of type integer.
|
||||
|
||||
Simard & L'Ecuyer (2011) [7].
|
||||
"""
|
||||
if np.isnan(n):
|
||||
return n # Keep the same type of nan
|
||||
if int(n) != n or n <= 0:
|
||||
return np.nan
|
||||
if x >= 1.0:
|
||||
return _select_and_clip_prob(1.0, 0.0, cdf=cdf)
|
||||
if x <= 0.0:
|
||||
return _select_and_clip_prob(0.0, 1.0, cdf=cdf)
|
||||
t = n * x
|
||||
if t <= 1.0: # Ruben-Gambino: 1/2n <= x <= 1/n
|
||||
if t <= 0.5:
|
||||
return _select_and_clip_prob(0.0, 1.0, cdf=cdf)
|
||||
if n <= 140:
|
||||
prob = np.prod(np.arange(1, n+1) * (1.0/n) * (2*t - 1))
|
||||
else:
|
||||
prob = np.exp(_log_nfactorial_div_n_pow_n(n) + n * np.log(2*t-1))
|
||||
return _select_and_clip_prob(prob, 1.0 - prob, cdf=cdf)
|
||||
if t >= n - 1: # Ruben-Gambino
|
||||
prob = 2 * (1.0 - x)**n
|
||||
return _select_and_clip_prob(1 - prob, prob, cdf=cdf)
|
||||
if x >= 0.5: # Exact: 2 * smirnov
|
||||
prob = 2 * scipy.special.smirnov(n, x)
|
||||
return _select_and_clip_prob(1.0 - prob, prob, cdf=cdf)
|
||||
|
||||
nxsquared = t * x
|
||||
if n <= 140:
|
||||
if nxsquared <= 0.754693:
|
||||
prob = _kolmogn_DMTW(n, x, cdf=True)
|
||||
return _select_and_clip_prob(prob, 1.0 - prob, cdf=cdf)
|
||||
if nxsquared <= 4:
|
||||
prob = _kolmogn_Pomeranz(n, x, cdf=True)
|
||||
return _select_and_clip_prob(prob, 1.0 - prob, cdf=cdf)
|
||||
# Now use Miller approximation of 2*smirnov
|
||||
prob = 2 * scipy.special.smirnov(n, x)
|
||||
return _select_and_clip_prob(1.0 - prob, prob, cdf=cdf)
|
||||
|
||||
# Split CDF and SF as they have different cutoffs on nxsquared.
|
||||
if not cdf:
|
||||
if nxsquared >= 370.0:
|
||||
return 0.0
|
||||
if nxsquared >= 2.2:
|
||||
prob = 2 * scipy.special.smirnov(n, x)
|
||||
return _clip_prob(prob)
|
||||
# Fall through and compute the SF as 1.0-CDF
|
||||
if nxsquared >= 18.0:
|
||||
cdfprob = 1.0
|
||||
elif n <= 100000 and n * x**1.5 <= 1.4:
|
||||
cdfprob = _kolmogn_DMTW(n, x, cdf=True)
|
||||
else:
|
||||
cdfprob = _kolmogn_PelzGood(n, x, cdf=True)
|
||||
return _select_and_clip_prob(cdfprob, 1.0 - cdfprob, cdf=cdf)
|
||||
|
||||
|
||||
def _kolmogn_p(n, x):
|
||||
"""Computes the PDF for the two-sided Kolmogorov-Smirnov statistic.
|
||||
|
||||
x must be of type float, n of type integer.
|
||||
"""
|
||||
if np.isnan(n):
|
||||
return n # Keep the same type of nan
|
||||
if int(n) != n or n <= 0:
|
||||
return np.nan
|
||||
if x >= 1.0 or x <= 0:
|
||||
return 0
|
||||
t = n * x
|
||||
if t <= 1.0:
|
||||
# Ruben-Gambino: n!/n^n * (2t-1)^n -> 2 n!/n^n * n^2 * (2t-1)^(n-1)
|
||||
if t <= 0.5:
|
||||
return 0.0
|
||||
if n <= 140:
|
||||
prd = np.prod(np.arange(1, n) * (1.0 / n) * (2 * t - 1))
|
||||
else:
|
||||
prd = np.exp(_log_nfactorial_div_n_pow_n(n) + (n-1) * np.log(2 * t - 1))
|
||||
return prd * 2 * n**2
|
||||
if t >= n - 1:
|
||||
# Ruben-Gambino : 1-2(1-x)**n -> 2n*(1-x)**(n-1)
|
||||
return 2 * (1.0 - x) ** (n-1) * n
|
||||
if x >= 0.5:
|
||||
return 2 * scipy.stats.ksone.pdf(x, n)
|
||||
|
||||
# Just take a small delta.
|
||||
# Ideally x +/- delta would stay within [i/n, (i+1)/n] for some integer a.
|
||||
# as the CDF is a piecewise degree n polynomial.
|
||||
# It has knots at 1/n, 2/n, ... (n-1)/n
|
||||
# and is not a C-infinity function at the knots
|
||||
delta = x / 2.0**16
|
||||
delta = min(delta, x - 1.0/n)
|
||||
delta = min(delta, 0.5 - x)
|
||||
|
||||
def _kk(_x):
|
||||
return kolmogn(n, _x)
|
||||
|
||||
return _derivative(_kk, x, dx=delta, order=5)
|
||||
|
||||
|
||||
def _kolmogni(n, p, q):
|
||||
"""Computes the PPF/ISF of kolmogn.
|
||||
|
||||
n of type integer, n>= 1
|
||||
p is the CDF, q the SF, p+q=1
|
||||
"""
|
||||
if np.isnan(n):
|
||||
return n # Keep the same type of nan
|
||||
if int(n) != n or n <= 0:
|
||||
return np.nan
|
||||
if p <= 0:
|
||||
return 1.0/n
|
||||
if q <= 0:
|
||||
return 1.0
|
||||
delta = np.exp((np.log(p) - scipy.special.loggamma(n+1))/n)
|
||||
if delta <= 1.0/n:
|
||||
return (delta + 1.0 / n) / 2
|
||||
x = -np.expm1(np.log(q/2.0)/n)
|
||||
if x >= 1 - 1.0/n:
|
||||
return x
|
||||
x1 = scu._kolmogci(p)/np.sqrt(n)
|
||||
x1 = min(x1, 1.0 - 1.0/n)
|
||||
|
||||
def _f(x):
|
||||
return _kolmogn(n, x) - p
|
||||
|
||||
return scipy.optimize.brentq(_f, 1.0/n, x1, xtol=1e-14)
|
||||
|
||||
|
||||
def kolmogn(n, x, cdf=True):
|
||||
"""Computes the CDF for the two-sided Kolmogorov-Smirnov distribution.
|
||||
|
||||
The two-sided Kolmogorov-Smirnov distribution has as its CDF Pr(D_n <= x),
|
||||
for a sample of size n drawn from a distribution with CDF F(t), where
|
||||
:math:`D_n &= sup_t |F_n(t) - F(t)|`, and
|
||||
:math:`F_n(t)` is the Empirical Cumulative Distribution Function of the sample.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
n : integer, array_like
|
||||
the number of samples
|
||||
x : float, array_like
|
||||
The K-S statistic, float between 0 and 1
|
||||
cdf : bool, optional
|
||||
whether to compute the CDF(default=true) or the SF.
|
||||
|
||||
Returns
|
||||
-------
|
||||
cdf : ndarray
|
||||
CDF (or SF it cdf is False) at the specified locations.
|
||||
|
||||
The return value has shape the result of numpy broadcasting n and x.
|
||||
"""
|
||||
it = np.nditer([n, x, cdf, None],
|
||||
op_dtypes=[None, np.float64, np.bool_, np.float64])
|
||||
for _n, _x, _cdf, z in it:
|
||||
if np.isnan(_n):
|
||||
z[...] = _n
|
||||
continue
|
||||
if int(_n) != _n:
|
||||
raise ValueError(f'n is not integral: {_n}')
|
||||
z[...] = _kolmogn(int(_n), _x, cdf=_cdf)
|
||||
result = it.operands[-1]
|
||||
return result
|
||||
|
||||
|
||||
def kolmognp(n, x):
|
||||
"""Computes the PDF for the two-sided Kolmogorov-Smirnov distribution.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
n : integer, array_like
|
||||
the number of samples
|
||||
x : float, array_like
|
||||
The K-S statistic, float between 0 and 1
|
||||
|
||||
Returns
|
||||
-------
|
||||
pdf : ndarray
|
||||
The PDF at the specified locations
|
||||
|
||||
The return value has shape the result of numpy broadcasting n and x.
|
||||
"""
|
||||
it = np.nditer([n, x, None])
|
||||
for _n, _x, z in it:
|
||||
if np.isnan(_n):
|
||||
z[...] = _n
|
||||
continue
|
||||
if int(_n) != _n:
|
||||
raise ValueError(f'n is not integral: {_n}')
|
||||
z[...] = _kolmogn_p(int(_n), _x)
|
||||
result = it.operands[-1]
|
||||
return result
|
||||
|
||||
|
||||
def kolmogni(n, q, cdf=True):
|
||||
"""Computes the PPF(or ISF) for the two-sided Kolmogorov-Smirnov distribution.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
n : integer, array_like
|
||||
the number of samples
|
||||
q : float, array_like
|
||||
Probabilities, float between 0 and 1
|
||||
cdf : bool, optional
|
||||
whether to compute the PPF(default=true) or the ISF.
|
||||
|
||||
Returns
|
||||
-------
|
||||
ppf : ndarray
|
||||
PPF (or ISF if cdf is False) at the specified locations
|
||||
|
||||
The return value has shape the result of numpy broadcasting n and x.
|
||||
"""
|
||||
it = np.nditer([n, q, cdf, None])
|
||||
for _n, _q, _cdf, z in it:
|
||||
if np.isnan(_n):
|
||||
z[...] = _n
|
||||
continue
|
||||
if int(_n) != _n:
|
||||
raise ValueError(f'n is not integral: {_n}')
|
||||
_pcdf, _psf = (_q, 1-_q) if _cdf else (1-_q, _q)
|
||||
z[...] = _kolmogni(int(_n), _pcdf, _psf)
|
||||
result = it.operands[-1]
|
||||
return result
|
||||
Some files were not shown because too many files have changed in this diff Show More
Reference in New Issue
Block a user