正态分布是很多計量資料比較分析的假設前提,是以在做比較分析之前要首先驗證樣本資料所代表的總體是否服從正态分布。當然對于比率資料的比較也需要滿足分布前提,通常是二項分布和泊松分布,對于二項分布的比率比較,一般不需要做分布的驗證。而對泊松分析的比率比較則需要事先驗證其分布,驗證方法就是卡方檢驗。
評價一進制正态性的方法有大緻分為以下幾種:
1、圖形驗證
- 直方圖
- 經驗累積機率圖(CDF)
- QQ圖
2、偏度和峰度
3、統計檢驗方法
常見的有以下幾種
- Kolmogorov-Smirnov檢驗
- Cramer-von Mises檢驗
- Anderson–Darling檢驗
- Shapiro–Wilk檢驗
一、圖形驗證
優點:直覺
缺點:1)、對于小樣本并不适用;2)、圖形方法以及下邊的峰度偏度法隻提供了一個粗糙而不正式的檢驗方法,沒有一個明确的決定準則。
1、頻率直方圖
直方圖需要注意的是資料的分區,對于小樣本資料分區過大或過小畫出來的直方圖可能都會影響效果,可以多試幾種分區,選擇其中較好的分區效果。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import norm
data=norm.rvs(loc=0,scale=1,size=1000) #生成1000個符合标準正态分布的浮點數。
#傳回的bins為清單,為頻率分區
n, bins, patches=plt.hist(data,
density=True, # true表示繪制的頻率直方圖總面積為1
bins=20) #預設分為10組,可以通過調節bins值,擴大分組。
#計算樣本的均值和标準差
mu=data.mean()
sigma=data.std()
#拟合一條最佳正态分布曲線y,y的值和分區bins的值對應
y = norm.pdf(bins, mu, sigma)
plt.plot(bins, y, 'r--') #繪制y的曲線
得到以下的頻率直方圖效果:
或者用seaborn繪制
import seaborn as sns
sns.set_palette("hls") #設定所有圖的顔色,使用hls色彩空間
sns.distplot(data,color="r",bins=20,kde=True)
plt.show()
輸出:
我們知道,标準正态分布的資料,大緻拟合一條“鐘形”的曲線。
2、經驗累積機率圖
第二種要介紹的圖是經驗累積機率圖,它的方法就是把數從小到大排列,求出每個值在總樣本中出現的機率,如樣本量是12,那麼每個數出現的機率就是1/12=0.083。所謂累積,就是從小到大開始算起,前i個數加起來出現的機率,如前10個數的累積機率是0.83,畫在圖中就是這樣的階梯形的線。
import numpy as np
import statsmodels.api as sm # recommended import according to the docs
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab
plt.rcParams['font.sans-serif'] = ['FangSong'] # 設定預設字型
plt.rcParams['axes.unicode_minus'] = False # 用于正常顯示'-'負号
#sample = np.random.uniform(0, 1, 50)
#非正态
#sample=np.array([1.7,2.8,1.4,2.6,-0.7,0.7,1.7,2.3,-28.00,2.4,0.0])
#正态
sample=np.array([9,18,13,0,5,-5,28,0,25,25,9,1])
ecdf = sm.distributions.ECDF(sample) #累計密度機率
#等差數列,用于繪制X軸資料
x = np.linspace(min(sample), max(sample)) #預設的個數為50個。
# x軸資料上值對應的累計密度機率
y = ecdf(x)
#繪制階梯圖
bins=plt.step(x, y,'b')
mu=sample.mean()
sigma=sample.std()
#建立一個包含1000個元素的樣本,提前設定均值和标準差
#loc表示均值
#scale表示标準差
#size表示個數
sample_normal=np.random.normal(loc=mu, scale=sigma, size=1000)
ecdf_normal = sm.distributions.ECDF(sample_normal) #累計密度機率
#等差數列,用于繪制X軸資料
x_normal = np.linspace(min(sample_normal), max(sample_normal),1000) #預設的個數為50個。
# x軸資料上值對應的累計密度機率
y_normal = ecdf_normal(x_normal)
#繪制階梯圖
plt.plot(x_normal,y_normal,'r--')
plt.title('樣本累計密度機率和正态累積機率密度')
plt.show()
圖中紅色的線是用樣本均值和标準差畫出來的理論曲線,它是用正态機率密度函數通過積分計算出來的。兩條曲線對比,如果挨得很近,則資料正态的可能性就很大。
注意這種圖不僅能夠直覺拟合正态曲線,也能拟合其它分布的曲線。
3、QQ圖
标準正态分布是一條“鐘形”的曲線;
分布同樣也是呈現“鐘形”,我們往往無法從頻率分布直方圖直接看出資料到底否服從正态分布呢?還是服從
分布?為了友善比較,我們建構三種資料集,分别服從正态分布、指數分布以及
分布,并繪制相應的直方圖,圖形如下:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import norm,t,lognorm,expon
plt.rcParams['font.sans-serif']=['SimHei'] #用來顯示中文标簽
plt.rcParams['axes.unicode_minus']=False #用來正常顯示負号
data={}
data['标準正态分布']=norm.rvs(loc=0,scale=1,size=1000) #生成1000個,滿足(0,1)的标準正态分布
data['指數分布']=expon.rvs(scale=10,size=1000) #生成1000個服從(1/λ,1/λ)的指數分布,scale為1/λ,指數分布的均值和标準差都等于1/λ。
data['學生t分布']=t.rvs(df=5,loc=0,scale=1,size=1000) #生成1000個服從t分布的資料,df表示學生t分布的自由度,詳細可查詢t分布和gamma分布的關系。
fig=plt.figure(figsize=(12,3))
for i,(name,data_) in enumerate(data.items()):
ax=fig.add_subplot(1,3,i+1)
n, bins, patches=plt.hist(data_,label=name,
density=True, # true表示繪制的頻率直方圖總面積為1
bins=20)
plt.title(name)
plt.show()
注意:
分布的機率密度函數為:
,其中
表示自由度,更多
分布和gamma分布的關系可見:《學生化t分布》,《統計學中,t分布中的自由度怎麼解釋?》。
資料分布可見下圖:
可以看出
分布比高斯(正态)分布有更長的尾巴(我們也叫其“厚尾”),也就是兩邊延伸得更開,這給出了
分布的一個重要性質——魯棒性,更長的尾巴意味着對于離群點能有更好的忍耐度,不會像高斯分布那樣敏感。
接下來我們來畫QQ圖,我們知道,對樣本資料标準化的公式為
,對其進行變形,
,若樣本資料近似于正态分布,在QQ圖上這些點近似的在直線
上,此直線的斜率是标準差
,截距為均值
。是以,利用QQ圖可以作直覺的正态性檢驗。若QQ圖上的點近似的在一條直線上,可以認為樣本資料來自正态總體。
- 繪制方法一:
#畫圖
plt.rcParams['font.sans-serif']=['SimHei'] #用來顯示中文标簽
plt.rcParams['axes.unicode_minus']=False #用來正常顯示負号
fig=plt.figure(figsize=(12,3))
x=np.array([-5,5])
for i,(name,data_) in enumerate(data.items()):
ax=fig.add_subplot(1,3,i+1)
mu=data_.mean() #樣本均值
sigma=data_.std() #樣本标準差
y=sigma*x+mu
plt.plot(x,y) #繪制截距為mu,斜率為sigma的直線
xs=norm.rvs(loc=0,scale=1,size=len(data_)) #抽取等樣本的标準正态值
xs.sort()
data_.sort()
plt.scatter(xs,data_,color='r')
plt.xlabel('标準正态值')
plt.ylabel('原始資料值')
plt.title(name)
plt.show()
輸出:
- 繪制方法二:
scipy中的stats有直接繪制QQ圖的方法,如下:
from scipy import stats
fig=plt.figure(figsize=(12,3))
for i,(name,data_) in enumerate(data.items()):
ax=fig.add_subplot(1,3,i+1)
stats.probplot(data_, dist=stats.norm, plot=ax)
plt.xlabel('标準正态值')
plt.ylabel('原始資料值')
plt.xlim(-4,4)
plt.title(name)
輸出:
根據QQ圖的形狀來判斷正态性
- 直線 正态
- 反“S”形 比正态厚尾
- “S”形 比正态薄尾
- 凸彎曲 右偏
- 凹彎曲 左偏
二、數值驗證
偏度講的是分布的偏斜程度,偏度<0,代表分布有左偏;偏度=0,表示分布左右對稱;偏度>0,表示分布有右偏,當然數值越大,偏斜程度也越大。
峰度講的是分布的尖銳度,或者說是胖還是瘦。峰度<0,則分布比較胖;峰度為>0,則分布比較瘦。
對于一組資料來說,如果計算出來的偏度和峰度都在0附近,那麼可以初步判斷其分布服從正态分布。
python有計算資料偏度和峰度的函數,如下,對上面的資料集求各自的偏度,峰度。
for i,(name,data_) in enumerate(data.items()):
data_s=pd.Series(data_)
print('%s偏度:%s;峰度:%s。'%(name,data_s.skew(),data_s.kurt()))
輸出:
标準正态分布偏度:-0.03483535284941679;峰度:0.008833727010514991。
指數分布偏度:1.930684125639356;峰度:4.741462640048853。
學生t分布偏度:0.845106666906179;峰度:5.391117437643347。
我們也可以調scipy中的stats()函數,計算不同分布下的偏度和峰度(輸出的均為總體的分布,不是樣本資料的分布)。
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"
norm.stats(loc=1,scale=2,moments='mvsk') #正态分布。
expon.stats(scale=10,moments='mvsk') #服從(1/λ,1/λ)的指數分布,scale為1/λ,指數分布的均值和标準差都等于1/λ。
t.stats(df=10000,loc=1,scale=2,moments='mvsk') #t分布,df代表自由度,df越大,越接近正态分布。
輸出:
(array(1.), array(4.), array(0.), array(0.))
(array(10.), array(100.), array(2.), array(6.))
(array(1.), array(4.00080016), array(0.), array(0.00060024))
四個參數分别代表均值、方差、偏度及峰度。
三、統計檢驗方法
以上驗證的方法都屬于直覺檢驗,樣本資料是否服從正态分布,還需要更加确切的檢驗。
常用的方法包括:
- Kolmogorov-Smirnov檢驗
- Cramer-von Mises檢驗
- Anderson–Darling檢驗
- Shapiro–Wilk檢驗
還有一些其它的檢驗方法,不一一列舉了。
1、Kolmogorov-Smirnov檢驗
正态性檢驗,最直覺的想法就是拿樣本資料與期望的理論分布進行對比,如果差異不大,則可以認為資料服從正态分布,Kolmogorov-Smirnov檢驗(下邊簡稱為“K-S檢驗”)方法就是這樣的。
K-S檢驗找出在每一個資料點上經驗累積機率與目标分布的累積機率之差的上界,列出公式是這樣的:
其中
代表樣本資料的經驗累積分布函數(cdf,上圖中的藍線),
代表與資料同均值、同方差的正态分布的累積分布函數(上圖中的紅線)。sup函數表示一組距離中的上确界,這是個數學概念,表示在原假設
的條件下,
的絕對值的最小上界。其意圖在于如果原假設成立,則
應該很小,如果很大,則原假設不成立。
但是,這個上确界怎麼求出來呢?請看下面的公式
其中k為樣本從小到大排列後的序數。從公式中看出
是經驗和目标累積機率之差和錯一位後再求出的差中最大的一個。
python實作
根據前面給的三個分布的資料,做K-S檢驗。
for i,(name,data_) in enumerate(data.items()):
test_stat =stats.kstest(data_,'norm')
print('%s:D統計量:%s;p_value值:%s。'%(name,test_stat[0],test_stat[1]))
輸出:
标準正态分布:D統計量:0.030836687076312264;p_value值:0.2925276653484524。
指數分布:D統計量:0.8027309998111607;p_value值:0.0。
學生t分布:D統計量:0.0494681189296266;p_value值:0.014470853698413954。
由結果我們知道,取α=0.05(顯著性水準),指數分布和學生t分布的的p_value值小于0.05,則拒絕原假設,即該分布不服從正态分布。
kstest 是一個很強大的檢驗子產品,除了正态性檢驗,還能檢驗 scipy.stats 中的其他資料分布類型,僅适用于連續分布的檢驗。
kstest(rvs, cdf, args=(), N=20, alternative='two_sided', mode='approx', **kwds)
對于正态性檢驗,我們隻需要手動設定三個參數即可:
- rvs:待檢驗的一組一維資料
- cdf:檢驗方法,例如'norm','expon','rayleigh','gamma',這裡我們設定為'norm',即正态性檢驗
- alternative:預設為雙尾檢驗,可以設定為'less'或'greater'作單尾檢驗
- model:'approx'(預設),使用檢驗統計量的精确分布的近視值,'asymp':使用檢驗統計量的漸進分布
注意:這種計算極差(理論和實際的最大差異)的方式,檢驗時對于樣本量的要求較高,适合于大樣本(一般
)資料的正态性檢驗。
2、Cramer-Von Mises檢驗(CVM檢驗)和Anderson–Darling檢驗(AD檢驗)
Cramer-Von Mises檢驗(CVM檢驗)
CVM檢驗的思想基本與K-S檢驗一緻,與K-S檢驗相比,CVM檢驗統計量度量的是經驗累積分布函數和目标累積分布函數的平方距離的積分。就是把每個資料點的差求平方以後相加,得到總的分布偏差,這樣就考慮了所有的差異點,而不是像K-S檢驗那樣隻考慮一個最大的。如下圖,其檢驗的統計量為:
其中
代表樣本資料的經驗累積分布函數(cdf,上圖中的藍線),
代表與資料同均值、同方差的正态分布的累積分布函數(上圖中的紅線)。檢驗的過程:在原假設
成立的條件下,
應該很小,如果
很大,則原假設不成立。
Anderson–Darling檢驗(AD檢驗)
AD檢驗實際上是對CVM檢驗的一種改進,其檢驗的統計量為:
AD檢驗實際上加上了一個權重函數
,這個權重,主要用于限制資料位置的權重,我們知道,如果資料落在左右兩邊的尾部位置,則
會特别小,無限就近了0,則
會變大,相當于在計算尾部資料的平方距離時,給予了更高的權重;計算中間資料的平方距離時,給予其較小的權重。結合下邊的圖,我們認為這樣設定是合理的。
3、Shapiro–Wilk檢驗(W檢驗)
K-S檢驗統計量用的是經驗累積機率與目标理論累積機率之差的最大值,這有點像計算極差;AD檢驗則把所有的差平方後求和,有點像計算方差。當然計算方差的穩健性(Robustness)要好一些,這可能也是AD檢驗得到更多應用的原因。接下來我們再來介紹一種檢驗方法——W檢驗。
W檢驗的統計量是:
表示第
個樣本次序統計量,是從小到大排序後的資料。
是标準正态分布中第
個樣本次序統計量标準化過後的值。即
。
我們可以将此代入到上面的公式得:
,此等式就類似我們求
和
的相關系數的公式。
可以把
看作是順序排列樣本值
和系數
之間相關系數的平方(squared correlation coefficient)或者是線性回歸的确定系數
(coefficient of determination
for linear regression),這個值最大為1,計算出的
值越接近1,樣本服從正态分布的可能性越大。當然也有判斷的臨界值,如果計算出的
值大于臨界值,則說明資料服從正态分布。
python實作
與 kstest 不同,shapiro 是專門用來做正态性檢驗的子產品。
for i,(name,data_) in enumerate(data.items()):
test_stat =stats.shapiro(data_)
print('%s:W統計量:%s;p_value值:%s。'%(name,test_stat[0],test_stat[1]))
輸出:
标準正态分布:W統計量:0.9982046484947205;p_value值:0.3791613280773163。
指數分布:W統計量:0.8033730387687683;p_value值:3.784428757365505e-33。
學生t分布:W統計量:0.9793980121612549;p_value值:1.0720994231272485e-10。
由結果我們知道,取α=0.05(顯著性水準),指數分布和學生
分布的的p_value值小于0.05,則拒絕原假設,即該分布不服從正态分布。
注意:
檢驗,精度較高,可适用于
的小樣本資料,有其他統計學家把其适用範圍擴充到5000,是以可以說
檢驗幾乎适用于所有的正态檢驗。正态性檢驗可優先考慮此方法。
了解更多正态性檢驗的推導,及如何選擇合适的正态性檢驗可以見《經典比較篇之三:正态分布檢驗方法(二)》;《經典比較篇之四:正态分布檢驗方法(三)》。
下邊截取其文章中部分内容,說明各種正态檢驗方法該如何選擇?
AD, RJ, 或KS: 哪一個正态檢驗最好?
前面介紹了幾種常用的正态分布檢驗方法,可能會有很多人感到困惑,哪種方法最好呢?應該如何選擇呢?理論上通常會進行功效上的對比,但實際應用上我們隻關注哪個相對更好一些,我們可以做各種模拟試驗來得出實際的結論,在這一點上已經有人做了,是以我就不費那個勁了,直接把别人的試驗成功翻譯出來,供大家參考。
Minitab軟體提供了三種正态檢驗:Anderson-Darling (AD), Ryan-Joiner (RJ), and Kolmogorov-Smirnov (KS)。AD檢驗是預設方法,但它是檢驗非正态的最佳方法嗎?讓我們比較一下在三種不同場景下這些檢驗正态的方法檢驗非正态資料的能力。我們用模拟的資料,但也基本反映了你在品質改進中分析資料時可能會碰到的常見情況。
場景1-制造過程會時不時會産生大的異常值,29個數模拟自均值=0,标準差=1的正态分布,1個數模拟自均值=0,标準差=4的正态分布。
場景2-制造過程發生了過程漂移,結果是分布上有平移。這産生了雙峰的分布,均值之間有4個σ的差距,樣本量為30。下圖展示了兩個正态分布之間均值存在的4個σ的差。
場景3-測量值實際服從非正态分布,典型的有失效時間和強度測量。本場景中,有30個資料模拟自Weibull分布(a=1,b=1.5)。
需要指出的是,三種場景并不是設計用來評估運用中心極限定理所做檢驗中正态假設的有效性,諸如1-樣本、2-樣本和配對t-檢驗。我們在這裡關注在用一個分布來估計制造缺陷的機率時檢驗出非正态的情形。
在場景1,Ryan-Joiner檢驗是個明顯的赢家,模拟結果見下圖。
【譯者注】縱坐标是拒絕正态的機率。
在場景2,Anderson-Darling檢驗最好,模拟結果如下。
在場景3,AD和RJ檢驗之間沒什麼不同。在檢驗非正态時兩個都比Kolmogorov-Smirnov檢驗更有效。模拟結果如下。
總結一下,Anderson-Darling檢驗總不是最壞的檢驗,但在檢驗一個4σ的異常值時不像RJ檢驗那樣有效。如果你要分析的制造過程資料有可能出現單個的異常值,Ryan-Joiner檢驗是最合适的。
RJ檢驗在兩種場景中表現很好,但當資料存在漂移時檢測非正态效果較差。如果你要分析的制造過程資料有可能由于未知的變化發生漂移,AD檢驗是最合适的。
KS檢驗在任何一種場景中都表現不好。
【譯者注】因為博文中沒有說明每個場景的模拟次數,是以模拟結論的可信性未知,另外這種模拟是否足夠嚴謹也是值得考慮的一個方面,大家參考着用。下面是作者的另一篇博文,關注的是資料舍入後正态性檢驗的變化。
正态檢驗和舍入
Jim Colton, 2013.12.4
所有的測量值都有舍入。在很多情況下,你不希望僅僅因為資料有舍入而拒絕正态。實際上,如果實際的分布是正态的,那麼正态分布是非常想要得到的資料模型,因為這樣可以把有舍入的測量值帶來的離散性弄得平滑一些。
有些正态檢驗在實際的分布服從正态時有非常大的比率因為舍入而拒絕(Anderson-Darling和Kolmogorov-Smirnov),而另一些似乎能忽略掉舍入(Ryan-Joiner和卡方)。
我們來考察一個極端的例子,資料如何在完美服從正态分布的情況下得到拒絕的結論,比如從一個均值為10和标準差為0.14的正态分布中抽一個樣本量為5000的樣本。
下面的圖展示了一部分舍入到小數點後兩位的資料及其一個直方圖和機率圖。直方圖和機率圖看起來非常好。Ryan-Joiner通過了正态檢驗,其p值大于0.10(左邊的機率圖)。然而Anderson-Darling檢驗的p值小于0.005(右邊的機率圖)。很明顯,在這個例子中拒絕正态是不合适的。
再模拟一個最常見的情況,n=30。模拟資料來自均值=0,标準差=1的正态分布,然後取整。下面就是這種模拟的一個例子。
在這種模拟的疊代過程中,Anderson-Darling的p值小于0.005,而Ryan-Joiner的p值大于0.10。
這種模拟的結果驚人地一緻,Anderson-Darling檢驗幾乎總是拒絕正态,Ryan-Joiner檢驗總是無法拒絕正态。模拟中也做了Kolmogorov-Smirnov和Chi-square檢驗。在避免因舍入而拒絕正态上,CS檢驗幾乎與RJ檢驗一樣的好。
第二個模拟考慮不那麼極端的舍入(舍入的程度取決于測量分辨力/标準差,比例越大,舍入越極端)。模拟資料來自均值=0和标準差=2的正态分布,然後舍入到最近的整數。其中的一個例子見下圖,在這種模拟的疊代中,Anderson-Darling的p值小于0.05,而Ryan-Joiner的p值大于0.10。
在第二個不太極端的舍入模拟中,AD和KS檢驗不那麼經常拒絕了。
第三個模拟與第二個模拟的舍入程度相同,但樣本量增大到n=100。由于樣本量較大,AD和KS檢驗又回到了幾乎總是拒絕的情況,RJ和CS檢驗也總是不因舍入而拒絕正态。
至此,我們隻讨論了因舍入而拒絕正态的情況,但RJ和CS能檢測到實際上是非正态的情況嗎?最後一個模拟與第一個模拟的舍入程度相同,但取自實際上非正态分布。模拟資料來自标準差為1的指數分布,然後舍入到最近的整數。AD和KS檢驗正确地拒絕正态,但可能是因為舍入,也可能是因為非正态。RJ和CS也能檢驗到非正态,但不像AD和KS那樣頻繁。
總結:RJ和CS避免因舍入(模拟1-3)而拒絕正态,而仍能檢測到來自完全非正态的資料(模拟4)。