1 檢查是否對(duì)稱
一般來(lái)說(shuō)夸赫,統(tǒng)計(jì)量較小的時(shí)候使用點(diǎn)圖混蔼,n較大的時(shí)候使用直方圖织阳,可以揭示一元分布的一個(gè)尾部比另一個(gè)長(zhǎng)的多的情況.
例子
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from scipy import stats
from matplotlib import style
mu, sigma = 0, 0.
s = np.random.normal(mu, sigma, 1000)
f, ax = plt.subplots(figsize=(10, 7))
plt.scatter(x=range(len(s)), y=s, c='r')
plt.xlim(0,500)
plt.show()
sns.distplot(s)
是不是很對(duì)稱啊
我們熟悉一下更多的正態(tài)分布樣圖
style.use('fivethirtyeight')
mu_params = [-1, 0, 1]
sd_params = [0.5, 1, 1.5]
x = np.linspace(-7, 7, 100)
f, ax = plt.subplots(len(mu_params), len(sd_params), sharex=True, sharey=True, figsize=(12,8))
for i in range(3):
for j in range(3):
mu = mu_params[i]
sd = sd_params[j]
y = stats.norm(mu, sd).pdf(x)
ax[i, j].plot(x, y)
ax[i, j].plot(0,0, label='mu={:3.2f}\nsigma={:3.2f}'.format(mu,sd), alpha=0)
ax[i, j].legend(fontsize=10)
ax[2,1].set_xlabel('x', fontsize=16)
ax[1,0].set_ylabel('pdf(x)', fontsize=16)
plt.suptitle('Gaussian PDF', fontsize=16)
plt.tight_layout()
plt.show()
2 區(qū)間檢查
一元正態(tài)分布屬于區(qū)間(μ-σ,μ+σ)內(nèi)的取值概率為0.683,
屬于區(qū)間(μ-2σ,μ+2σ)內(nèi)的取值概率為0.954
在(μ-3σ,μ+3σ)內(nèi)的取值概率為0.997
sigma = 2
list(map(lambda x: stats.norm.cdf(x*sigma,loc=0,scale=sigma) - stats.norm.cdf(-x*sigma,loc=0,scale=sigma), range(1, 6)))
[0.6826894921370859,
0.9544997361036416,
0.9973002039367398,
0.9999366575163338,
0.9999994266968562]
3 QQ圖與PP圖
QQ圖(Quantile Quantile Plot)有兩個(gè)作用秆乳,1檢查一組數(shù)據(jù)是否服從同一分布般妙,2檢查兩個(gè)分布是否服從同一分布
QQ圖原理是比較兩組數(shù)據(jù)的累計(jì)分布函數(shù)來(lái)判斷兩組數(shù)據(jù)是否是服從同一分布纪铺,所以第一步我們應(yīng)該做兩組數(shù)據(jù)的累計(jì)分布相速。首先碟渺,作為對(duì)比我們看一下標(biāo)準(zhǔn)正太分布的累計(jì)分布圖
x = np.linspace(-5, 5, 100)
y = stats.norm.cdf(x, 0, 1)
plt.plot(x, y)
然后,繪制目標(biāo)數(shù)據(jù)(這里使用隨機(jī)生成的數(shù)據(jù)集)的累計(jì)分布函數(shù)圖
x = np.random.normal(0, 1, 100)
sorted_ = np.sort(x)
yvals = np.arange(len(sorted_))/float(len(sorted_))
plt.plot(sorted_, yvals)
直觀上對(duì)比突诬,目標(biāo)累計(jì)分布函數(shù)圖和標(biāo)準(zhǔn)正太累計(jì)分布函數(shù)圖差異不大苫拍,事實(shí)是不是這樣呢?最后我們就可以做pp圖做對(duì)比旺隙。
x_label = stats.norm.ppf(yvals) #對(duì)目標(biāo)累計(jì)分布函數(shù)值求標(biāo)準(zhǔn)正太分布累計(jì)分布函數(shù)的逆
plt.scatter(x_label, sorted_)
stats.probplot(x, dist="norm", plot=plt)
plt.show()
下面做個(gè)QQ圖進(jìn)行比較
import statsmodels.api as sm
nsample = 100
plt.figure(figsize=(10, 8))
sm.qqplot(stats.t.rvs(25, size=nsample), line='45')
除了正態(tài)分布绒极,我們還可以檢測(cè)t分布(自動(dòng)度較小)
nsample = 100
plt.figure(figsize=(10, 8))
x = stats.t.rvs(3, size=nsample)
res = stats.probplot(x, plot=plt)
檢測(cè)t分布(自動(dòng)度較大)
nsample = 100
plt.figure(figsize=(10, 8))
x = stats.t.rvs(25, size=nsample)
res = stats.probplot(x, plot=plt)
兩種正態(tài)分布的混合
nsample = 100
plt.figure(figsize=(10, 8))
x = stats.norm.rvs(loc=[0,5], scale=[1,1.5],size=(nsample,2)).ravel()
res = stats.probplot(x, plot=plt)
loggamma分布
plt.figure(figsize=(10, 8))
x = stats.loggamma.rvs(c=2.5, size=500)
stats.probplot(x, dist=stats.loggamma, sparams=(2.5,), plot=plt)
4 正態(tài)的相關(guān)系數(shù)
例如序列x = [-1.0, -0.1, 0.16, 0.41, 0.62, 0.8, 1.26, 1.54, 1.71, 2.3]
則均值sum(x)/len(x) = 0.77
概率水平
n =10
list(map(lambda x:(x-0.5)/n, range(1, n+1)))
[0.05, 0.15, 0.25, 0.35, 0.45, 0.55, 0.65, 0.75, 0.85, 0.95]
對(duì)應(yīng)的分位數(shù)為
map(lambda x:stats.norm.ppf((x-0.5)/n), range(1, n+1))
為了求相關(guān)系數(shù),先求三個(gè)數(shù),然后計(jì)算rq
x = [-1.0, -0.1, 0.16, 0.41, 0.62, 0.8, 1.26, 1.54, 1.71, 2.3]
n=len(x)
x_ = sum(x)/n
r1 = sum(map(lambda x: (x[0]-x_)*x[1], zip(x, map(lambda x:stats.norm.ppf((x-0.5)/n), range(1, n+1)))))
r2 = sum(map(lambda x: (x-x_)*(x-x_), x))
r3 = sum(map(lambda x:stats.norm.ppf((x-0.5)/n)*stats.norm.ppf((x-0.5)/n), range(1, n+1)))
rq = r1/(np.sqrt(r2)*np.sqrt(r3))
rq = 0.9943587408592187
在顯著性水平10%下蔬捷,把rq = 0.9943587408592187與查表n=10, a=0.10相對(duì)照垄提,進(jìn)行正態(tài)性檢驗(yàn),則rq>0.9341 周拐,因此不拒絕正態(tài)性假定
5 柯?tīng)柲陕宸?斯米洛夫檢驗(yàn)Kolmogorov-Smirnov test
柯?tīng)柲缏宸?斯米爾諾夫檢驗(yàn)(以下簡(jiǎn)稱K-S檢驗(yàn))是用累計(jì)次數(shù)或累計(jì)頻率來(lái)判斷兩組數(shù)據(jù)之間是否存在顯著差異的方法铡俐。它是將需要做統(tǒng)計(jì)分析的數(shù)據(jù)和另一組標(biāo)準(zhǔn)數(shù)據(jù)進(jìn)行對(duì)比,求得它和標(biāo)準(zhǔn)數(shù)據(jù)之間的偏差的方法妥粟。
藍(lán)線表示數(shù)據(jù)审丘,紅線表示假設(shè)假定符合的分布。
而X軸表示數(shù)據(jù)值的大小勾给,Y軸表示的數(shù)據(jù)累計(jì)所占百分比滩报。如果簡(jiǎn)單理解實(shí)際上就是概率密度函數(shù)的積分。這這個(gè)圖里面紅線實(shí)際上就是正態(tài)分布的情況播急,而藍(lán)線因?yàn)槭请x散化的數(shù)據(jù)脓钾,所以呈現(xiàn)的是階梯狀。兩條線之間的最大距離也就是黑色箭頭表現(xiàn)的位置就表示了二者之間的最大區(qū)別程度桩警,稱為D可训,D值的大小則決定了兩組數(shù)據(jù)間的差異。用這種百分?jǐn)?shù)的差別來(lái)表現(xiàn)差異,有一個(gè)最明顯的好處沉噩,那就是不會(huì)因?yàn)槟骋粋€(gè)點(diǎn)的異常而否定所以的點(diǎn)捺宗。此外,KEST還可以檢驗(yàn)多種分布川蒙,只需要把紅線換成其他的線即可蚜厉。
到這里我們僅僅得到了D值,還不能完全判定兩者的符合程度畜眨,這時(shí)候還需要引入顯著度α(alpha默認(rèn)為0.05)
kstest(x, 'norm')
KstestResult(statistic=0.36355946289143287, pvalue=0.1084952486818282)
由于pvalue=0.1084952486818282 > 0.05 因此不能拒絕原來(lái)假設(shè)
6 夏皮羅-威爾克檢驗(yàn)
該檢驗(yàn)是由S.S.Shapiro與M.B.Wilk提出的昼牛,又被稱之為W檢驗(yàn),主要檢驗(yàn)研究對(duì)象是否符合正態(tài)分布康聂。假設(shè): 一定樣本量n(8<n<50)的研究對(duì)象總是符合正態(tài)分布贰健。
將樣本量為n的樣本按照大小順序編排,然后根據(jù)公式計(jì)算統(tǒng)計(jì)量W的值恬汁,該值越接近于1伶椿,且顯著水平大于0.05時(shí),我們就沒(méi)法拒絕原假設(shè)
xi為排序后的樣本數(shù)據(jù),ai為待估常量,統(tǒng)計(jì)量越大則表示數(shù)據(jù)越符合正態(tài)分布氓侧,但是僅憑這一個(gè)參數(shù)是不夠的脊另,在非正態(tài)分布的小樣本數(shù)據(jù)中也經(jīng)常會(huì)出現(xiàn)較大的W值。該統(tǒng)計(jì)量的分布是未知的约巷,因此需要通過(guò)模擬或者查表來(lái)估計(jì)其概率偎痛。由于原假設(shè)是其符合正態(tài)分布,所以當(dāng)P值小于指定顯著水平時(shí)表示其不符合正態(tài)分布独郎。
x = [-1.0, -0.1, 0.16, 0.41, 0.62, 0.8, 1.26, 1.54, 1.71, 2.3]
stats.shapiro(x)
(0.9899559020996094, 0.9967610239982605)
0.9967610239982605>0.05 妥妥的正態(tài)分布
7 正態(tài)性檢測(cè)匯總
from statsmodels.stats.diagnostic import lillifors
group1=[2,3,7,2,6]
group2=[10,8,7,5,10]
group3=[10,13,14,13,15]
list_groups=[group1,group2,group3]
list_total=group1+group2+group3
#正態(tài)分布測(cè)試
def check_normality(testData):
#20<樣本數(shù)<50用normal test算法檢驗(yàn)正態(tài)分布性
if 20<len(testData) <50:
p_value= stats.normaltest(testData)[1]
if p_value<0.05:
print("use normaltest")
print("data are not normal distributed")
return False
else:
print("use normaltest")
print("data are normal distributed")
return True
#樣本數(shù)小于50用Shapiro-Wilk算法檢驗(yàn)正態(tài)分布性
if len(testData) <50:
p_value= stats.shapiro(testData)[1]
if p_value<0.05:
print("use shapiro:")
print("data are not normal distributed")
return False
else:
print("use shapiro:")
print("data are normal distributed")
return True
if 300>=len(testData) >=50:
p_value= lillifors(testData)[1]
if p_value<0.05:
print("use lillifors:")
print("data are not normal distributed")
return False
else:
print("use lillifors:")
print("data are normal distributed")
return True
if len(testData) >300:
p_value= stats.kstest(testData,'norm')[1]
if p_value<0.05:
print("use kstest:")
print("data are not normal distributed")
return False
else:
print( "use kstest:")
print("data are normal distributed")
return True
#對(duì)所有樣本組進(jìn)行正態(tài)性檢驗(yàn)
def NormalTest(list_groups):
for group in list_groups:
#正態(tài)性檢驗(yàn)
status=check_normality(group1)
if status==False :
return False
#對(duì)所有樣本組進(jìn)行正態(tài)性檢驗(yàn)
NormalTest(list_groups)