SciPy - 正态性 与 KS 检验 与 boxcox正态化

假设检验实质上是对原假设是否正确进行检验，因此检验过程中要使原假设得到维护，使之不轻易被拒绝；否定原假设必须有充分的理由。同时，当原假设被接受时，也只能认为否定该假设的根据不充分，而不是认为它绝对正确

def kstest(rvs, cdf, args=(), N=20, alternative='two-sided', mode='approx'):
    """
    Parameters
    ----------
    rvs : str, array_like, or callable
        If a string, it should be the name of a distribution in `scipy.stats`.
        If an array, it should be a 1-D array of observations of random
        variables.
        If a callable, it should be a function to generate random variables;
        it is required to have a keyword argument `size`.
    Returns
    -------
    statistic : float
        KS test statistic, either D, D+ or D-.
    pvalue :  float
        One-tailed or two-tailed p-value.
    """

######### 非正态 #########
x1 = np.linspace(-15, 15, 9)
print(kstest(x1, 'norm'))
# KstestResult(statistic=0.4443560271592436, pvalue=0.03885014270517116)  <0.05

######### 正态 #########
np.random.seed(1000)
x2 = np.random.randn(100)
print(kstest(x2, 'norm'))
# KstestResult(statistic=0.06029093862878099, pvalue=0.8604070909241421)   >0.05

#### 同一个数据(服从正态分布)，不同的参数有截然不同的检测结果，说明 ks 检测 正态性 比较麻烦
x3 = np.random.normal(100, 0.01, 1000)
print(kstest(x3, 'norm'))
# KstestResult(statistic=1.0, pvalue=0.0)           <0.05
print(kstest(x3, 'norm', alternative='greater'))
# KstestResult(statistic=0.0, pvalue=1.0)           >0.05

import numpy as np
from scipy.stats import ks_2samp, kstest

beta=np.random.beta(7, 5, 1000)
norm=np.random.normal(0, 1, 1000)
print(ks_2samp(beta, norm))
# Ks_2sampResult(statistic=0.578, pvalue=7.844864182954565e-155)

# p-value比指定的显著水平（假设为5%）小，则我们完全可以拒绝假设：beta和norm不服从同一分布

def shapiro(x):
    """
    Parameters
    ----------
    x : array_like
        Array of sample data.

    Returns
    -------
    W : float
        The test statistic.
    p-value : float
        The p-value for the hypothesis test.
    """

######### 非正态 #########
x = np.random.rand(100)
print(stats.shapiro(x))
# (0.9387611746788025, 0.00016213695926126093)   <0.05

######### 正态 #########
x1 = stats.norm.rvs(loc=5, scale=3, size=100)
print(stats.shapiro(x1))
# (0.9818177223205566, 0.18372832238674164)     >0.05

np.random.seed(1000)
x2 = np.random.randn(100)
print(stats.shapiro(x2))
# (0.9930729269981384, 0.89242023229599)        >0.05

#### 同样的数据，ks 检测 非正态，shapiro 检测 正态
x3 = np.random.normal(100, 0.01, 1000)
print(stats.shapiro(x3))
# (0.9976787567138672, 0.17260634899139404)     >0.05

def normaltest(a, axis=0, nan_policy='propagate'):
    """
    Parameters
    ----------
    a : array_like
        The array containing the sample to be tested.
    Returns
    -------
    statistic : float or array
    pvalue : float or array
       A 2-sided chi squared probability for the hypothesis test.
    """

######### 非正态 #########
x = np.random.rand(100)
print(normaltest(x))
# NormaltestResult(statistic=17.409796250892015, pvalue=0.00016577184786199797) <0.05

######### 正态 #########
x1 = np.random.randn(10, 20)
print(normaltest(x1, axis=None))
#NormaltestResult(statistic=1.8245865103063612, pvalue=0.4016021909152733)  >0.05

np.random.seed(1000)
x2 = np.random.randn(100)
print(normaltest(x2))
# NormaltestResult(statistic=1.2417269613653144, pvalue=0.5374801334521462)       >0.05


#### 同样的数据，ks 检测 非正态，shapiro 检测 正态, normaltest 检测 正态
x3 = np.random.normal(100, 0.01, 1000)
print(normaltest(x3))
# NormaltestResult(statistic=3.3633106484154944, pvalue=0.18606572188584633)    >0.05

def anderson(x, dist='norm'):
    """
    Anderson-Darling test for data coming from a particular distribution.
    x : array_like
        Array of sample data.
    dist : {'norm','expon','logistic','gumbel','gumbel_l', gumbel_r',
        'extreme1'}, optional
        the type of distribution to test against.  The default is 'norm'
        and 'extreme1', 'gumbel_l' and 'gumbel' are synonyms.

    Returns
    -------
    statistic : float
        The Anderson-Darling test statistic.
    critical_values : list
        The critical values for this distribution.
    significance_level : list
        The significance levels for the corresponding critical values
        in percents.  The function returns critical values for a
        differing set of significance levels depending on the
        distribution that is being tested against.
    """

# If the returned statistic is larger than these critical values then for the corresponding significance level,
# the null hypothesis that the data come from the chosen distribution can be rejected.

######### 非正态 #########
x1 = np.random.rand(100)
print(anderson(x1))
# AndersonResult(statistic=2.297910919361925, critical_values=array([0.555, 0.632, 0.759, 0.885, 1.053]), significance_level=array([15. , 10. ,  5. ,  2.5,  1. ]))
### 统计数 大于 评判值，拒绝假设，非正态

######### 正态 #########
x2 = np.random.randn(100)
print(anderson(x2))
# AndersonResult(statistic=0.35148092757619054, critical_values=array([0.555, 0.632, 0.759, 0.885, 1.053]), significance_level=array([15. , 10. ,  5. ,  2.5,  1. ]))
### 统计数 小于 评判值，正态

#### 同样的数据，ks 检测 非正态，shapiro 检测 正态
np.random.seed(1000)
x3 = np.random.normal(100, 0.01, 1000)
print(anderson(x3))
# AndersonResult(statistic=0.201867508676969, critical_values=array([0.574, 0.653, 0.784, 0.914, 1.088]), significance_level=array([15. , 10. ,  5. ,  2.5,  1. ]))
### 统计数 小于 评判值，正态

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as st

fig, axs = plt.subplots(3, 1)

######### 非正态 #########
x1 = np.random.randint(1, 100, 100)
st.probplot(x1, plot=axs[0])

######### 正态 #########
np.random.seed(1000)
x2 = np.random.randn(100)
st.probplot(x2, plot=axs[1])

#### ks 检验非正态，看看 qq 图
x3 = np.random.normal(100, 0.01, 1000)
st.probplot(x3, plot=axs[2])

plt.show()

from scipy.stats import kstest, shapiro
from scipy.stats import boxcox
import matplotlib.pyplot as plt
from scipy import stats

fig, axs = plt.subplots(2, 1)
np.random.seed(12345)

x = np.random.normal(100, 20, 200) + np.random.normal(40, 5, 200)
x = np.log2(x)

print(shapiro(x))
# (0.9857848882675171, 0.04184681549668312)     ### 非正态
stats.probplot(x, plot=axs[0])

x2, _ = boxcox(x)
print(shapiro(x2))
# (0.9948143362998962, 0.7225888967514038)      ### 正态
stats.probplot(x2, plot=axs[1])
plt.show()