Python統(tǒng)計學包scipystats手冊

1、Statistics(scipy.stats)Statistics(scipy.stats)1介紹2隨機變量2獲得幫助2通用方法4位移與縮放6形態(tài)參數(shù)7凍結(jié)分布9廣播9離散分布的特殊之處10分布擬合12性能問題與注意事項12遺留問題12構造具體的分布13創(chuàng)建一個連續(xù)分布，繼承rv_continuous類13繼承rv_discrete類14樣本分析19描述統(tǒng)計20T檢驗和KS檢驗21分布尾部23正態(tài)分布的特殊檢驗25比較兩個樣本27均值27對于兩個不同的樣本進行的KS檢僉28核密度估計28單元估計28多元估計36介紹c:iknowdocsharedatacurworkreferencetutor

2、ialstats.html在這個教程我們討論一部分scipy.stats模塊的特性。這里我們的意圖是提供給使用者一個關于這個包的實用性知識。我們推薦referencemanual來介紹更多的細節(jié)。注意：這個文檔還在發(fā)展中。隨機變量有些通用的概率分布類被圭寸裝在continuousrandomvariables以及discreterandomvariables中。有80多個連續(xù)性隨機變量(RVs)以及10余個離散隨機變量已經(jīng)用這些類建立。同樣，新的程序和分布可以被用戶新建(如果你構造了一個，請?zhí)峁┧o我們幫助發(fā)展這個包)所有統(tǒng)計函數(shù)被放在子包scipy.stats中，且有這些函數(shù)的一個幾乎完整的

3、列表可以使用info(stats)獲得。這個列表里的隨機變量也可以從stats子包的docstring中獲得介紹。在接下來的討論中，我們著重于連續(xù)性隨機變量(RVs)o幾乎所有離散變量也符合下面的討論，但是我們也要指出一些區(qū)別在“離散分布的特殊之處”中。獲得幫助所有分布可以使用help函數(shù)得到解釋。為獲得這些信息只需要使用像這樣的簡單調(diào)用:>>>>>>fromscipy import stats>>>fromscipy.stats import norm>>>printnorm . doc作為例子，我們用這種方式找分布的上

4、下界>>>>>>print'boundsofdistributionlower:%s,upper:%s'%(norm.a,norm.b)boundsofdistributionlower:-inf,upper:inf我們可以通過調(diào)用dir(norm)來獲得關于這個(正態(tài))分布的所有方法和屬性。應該看到，一些方法是私有方法盡管其并沒有以名稱表示出來(比如它們前面沒有以下劃線開頭)，比如veccdf就只用于內(nèi)部計算(試圖使用那些方法將引發(fā)警告，因為它們可能會在后續(xù)開發(fā)中被移除)為了獲得真正的主要方法，我們列舉凍結(jié)分布的方法(我們將在下文解釋何謂“

5、凍結(jié)分布”)>>>>>>rv=norm()>>>dir(rv)#reformatted'class','delattr','dict','doc','getattribute','hash','init','module','new','reduce','reduceex','repr_','_setattr_','_str_&

6、#39;,'_weakref','args','cdf','dist','entropy','isf','kwds','moment','pdf','pmf','ppf','rvs','sf','stats'最后，我們能通過內(nèi)省獲得所有的可用分布的信息。>>>>>>importwarnings>>>warnings.si

7、mplefilter('ignore',DeprecationWarning)>>>dist_continu=dfordindir(stats)if.isinstance(getattr(stats,d),stats.rv_continuous)>>>dist_discrete=dfordindir(stats)if.isinstance(getattr(stats,d),stats.rv_discrete)>>>print'numberofcontinuousdistributions:',len(dist

8、_continu)numberofcontinuousdistributions:84>>>print'numberofdiscretedistributions:',len(dist_discrete)numberofdiscretedistributions:12通用方法連續(xù)隨機變量的主要公共方法如下：?rvs:隨機變量(就是從這個分布中抽一些樣本)?pdf:概率密度函數(shù)。?cdf:累計分布函數(shù)?sf:殘存函數(shù)(1-CDF)?ppf:分位點函數(shù)(CDF的逆)?isf:逆殘存函數(shù)(sf的逆)?stats:返回均值，方差，(費舍爾)偏態(tài)，(費舍爾)峰度。?mo

9、ment:分布的非中心矩。讓我們使用一個標準的RV作為例子。>>>>>>norm.cdf(0)0.5為了計算在一個點上的cdf,我們可以傳遞一個列表或一個numpy數(shù)組>>>> >>norm.cdf(-1.,0,1)array(0.15865525,0.5,0.84134475)> >>importnumpyasnp> >>norm.cdf(np.array(-1.,0,1)array(0.15865525,0.5,0.84134475)相應的，像pdf,cdf之類的簡單方法可以用np.

10、vectorize矢量化.其他實用的方法可以像這樣使用。>>>> >>norm.mean(),norm.std(),norm.var()(0.0,1.0,1.0)> >>norm.stats(moments="mv")(array(0.0),array(1.0)為了找到一個分布的中心，我們可以使用分位數(shù)函數(shù)ppf,它是cdf的逆>>>> >>norm.ppf(0.5)0.0為了產(chǎn)生一個隨機變量集合。>>>>>>norm .rvs(size = 5)a

11、rray(-0.35687759,1.34347647,-0.11710531,-1.00725181,-0.51275702)不要認為norm.rvs(5)產(chǎn)生了五個變量。>>>> >>norm.rvs(5)7.5814欲知其意，請看下一部分的內(nèi)容。位移與縮放所有連續(xù)分布可以操縱loc以及scale參數(shù)作為修正location和scale的方式。作為例子，標準正態(tài)分布的location是均值而scale是標準差。>>>> >>norm.stats(loc=3,scale=4,moments="mv")

12、(array(3.0),array(16.0)通常經(jīng)標準化的分布的隨機變量X可以通過變換(X-loc)/scale獲得。它們的默認值是loc=0以及scale=1.聰明的使用loc與scale可以幫助以靈活的方式調(diào)整標準分布達到所想要的效果。為了進一步說明縮放的效果，下面給出期望為1/人指數(shù)分布的cdfoF(x)=1-exp(-入x)通過像上面那樣使用scale,可以看到得到想要的期望值。>>>> >>fromscipy.statsimportexpon> >>expon.mean(scale=3.)3.0均勻分布也是令人感興趣的：>

13、>>> >>fromscipy.statsimportuniform> >>uniform.cdf(0,1,2,3,4,5,loc=1,scale=4)array(0.,0.,0.25,0.5,0.75,1.)最后，聯(lián)系起我們在前面段落中留下的norm.rvs(5)的問題。事實上，像這樣調(diào)用一個分布，其第一個參數(shù)，像之前的5,是把loc參數(shù)調(diào)到了5,讓我們看：>>>> >>np.mean(norm.rvs(5,size=500)4.983550784784704在這里，為解釋最后一段的輸出：norm.rvs(5

14、)產(chǎn)生了一個正態(tài)分布變量，其期望，即loc=5.我傾向于明確的使用loc,scale作為關鍵字而非像上面那樣利用參數(shù)的順序。因為這看上去有點令人困惑。我們在我們解釋“凍結(jié)RV”的主題之前澄清這一點。形態(tài)參數(shù)雖然一個一般的連續(xù)隨機變量可以通過賦予loc和scale參數(shù)確定，但一些分布還需要額外的形態(tài)參數(shù)。作為例子，看到這個伽馬分布，這是它的密度函數(shù)Y(x,a)=入(入x)a-ir(a)e-入x,要求一個形態(tài)參數(shù)a。注意到人的設置可以通過設置scale關鍵字為1/人進行。讓我們檢查伽馬分布的形態(tài)參數(shù)的名字的數(shù)量。(我們知道從上面知道其應該為1)>>>>>>fro

15、mscipy.statsimportgamma>>>gamma.numargs1>>>gamma.shapes'a'現(xiàn)在我們設置形態(tài)變量的值為1以變成指數(shù)分布。所以我們可以容易的比較是否得到了我們所期望的結(jié)果。>>>>>>gamma(1,scale=2.).stats(moments="mv")(array(2.0),array(4.0)注意我們也可以以關鍵字的方式指定形態(tài)參數(shù)：>>>>>>gamma(a=1,scale=2.).stats(momen

16、ts="mv")(array(2.0),array(4.0)凍結(jié)分布不斷地傳遞loc與scale關鍵字最終會讓人厭煩。而凍結(jié)RV的概念被用來解決這個問題。>>>>>>rv=gamma(1,scale=2.)通過使用rv,在任何情況下我們不再需要包含scale與形態(tài)參數(shù)。顯然，分布可以被多種方式使用，我們可以通過傳遞所有分布參數(shù)給對方法的每次調(diào)用(像我們之前做的那樣)或者可以對一個分布對象凍結(jié)參數(shù)。讓我們看看是怎么回事：>>>>>>rv.mean(),rv.std()(2.0,2.0)這正是我們應該得到的

17、。廣播c:iknowdocsharedatacurworkipyreferencetutorialstats.html像pdf這樣的簡單方法滿足numpy的廣播規(guī)則。作為例子，我們可以計算t分布的右尾分布的臨界值對于不同的概率值以及自由度。>>>> >>stats.t.isf(0.1,0.05,0.01,10,11)array(1.37218364,1.81246112,2.76376946,1.36343032,1.79588482,2.71807918)這里，第一行是以10自由度的臨界值，而第二行是以11為自由度的臨界值。所以，廣播規(guī)則與下面調(diào)用了兩次i

18、sf產(chǎn)生的結(jié)果相同>>>> >>stats.t.isf(0.1,0.05,0.01,10)array(1.37218364,1.81246112,2.76376946)> >>stats.t.isf(0.1,0.05,0.01,11)array(1.36343032,1.79588482,2.71807918)但是如果概率數(shù)組，如0.1,0.05,0.01與自由度數(shù)組，如10,11,12具有相同的數(shù)組形態(tài)，則進行對應匹配，我們可以分別得到10%,5%,1%尾的臨界值對于10,11,12的自由度。>>>> >&g

19、t;stats.t.isf(0.1,0.05,0.01,10,11,12)array(1.37218364,1.79588482,2.68099799)離散分布的特殊之處c:iknowdocsharedatacurworkorialstats.html離散分布的方法的大多數(shù)與連續(xù)分布很類似。當然像pdf被更換為密度函數(shù)pmf,沒有估計方法，像fit就不能用了。而scale不是一個合法的關鍵字參數(shù)。Location參數(shù)，關鍵字loc則仍然可以使用用于位移。cdf的計算要求一些額外的關注。在連續(xù)分布的情況下，累積分布函數(shù)在大多數(shù)標準情況下是嚴格遞增的，所以有唯一的逆。而cdf在離散分布，無論如何，

20、是階躍函數(shù)，所以cdf的逆，分位點函數(shù)，要求一個不同的定義：ppf(q)=minx:cdf(x)>=q,xinteger為了更多信息可以看這里。我們可以看這個超幾何分布的例子>>>> >>fromscipy.statsimporthypergeom> >>M,n,N=20,7,12如果我們在一些整數(shù)點使用cdf,則它們的cdf值再作用ppf會回到開始的值>>>> >>x=np.arange(4)*2>>>xarray(0,2,4,6)> >>prb=hyperge

21、om.cdf(x,M,n,N)>>>prbarray(0.44066,0.05211,0.69301,0.9897832817337386)> >>hypergeom.ppf(prb,M,n,N)array(0.,2.,4.,6.)如果我們使用的值不是cdf的函數(shù)值，則我們得到一個更高的值。>>>> >>hypergeom.ppf(prb+1e-8,M,n,N)array(1.,3.,5.,7.)>>>hypergeom.ppf(prb-1e-8,M,n,N)array(0.,2.,4.,6.)分布擬合非

22、凍結(jié)分布的參數(shù)估計的主要方法：fit:分布參數(shù)的極大似然估計，包括location與scalefit_loc_scale:給定形態(tài)參數(shù)確定下的location和scale參數(shù)的估計nnlf:負對數(shù)似然函數(shù)expect:計算函數(shù)pdf或pmf的期望值。性能問題與注意事項每個方法的性能與運行速度根據(jù)分布的不同表現(xiàn)差異極大。方法的結(jié)果可以由兩種方式獲得，精確計算或使用獨立于各具體分布的通用算法。對于精確計算，一個分布能使用這種方式的第一種情況，這個分布是包中直接給你的（如標準正態(tài)分布），第二，給出解析形式，第三通過scipy.special或numpy.special或numpy.random的rv

23、s特殊函數(shù)給出。一般來說使用精確計算會比較快。另一方面，通用方法被用于當分布沒有被指派明確的計算方法時使用。為了定義一個分布，只有pdf或cdf是必須的；通用方法使用數(shù)值積分和求根法得到結(jié)果。作為例子，rgh=stats.gausshyper.rvs（0.5,2,2,2,size=100）以這種方式創(chuàng)建了100個隨機變量（抽了100個值），這在我的電腦上花了19秒（譯者：我花了3.5秒），對比取一百萬個標準正態(tài)分布的值只需要1秒。遺留問題scipy.stats里的分布最近進行了升級并且被仔細的檢查過了，不過仍有一些問題存在。?分布在很多參數(shù)區(qū)間上的值被測試過了，然而在一些奇葩的邊界，仍然可能有

24、錯誤的值存在。?fit的極大似然估計以默認值作為初始參數(shù)將不會工作的很好，用戶必須指派合適的初始參數(shù)。并且，對于一些分布使用極大似然估計本身就不是一個好的選擇。構造具體的分布下一個例子展示了如何建立你自己的分布。更多的例子見分布用法以及統(tǒng)計檢驗創(chuàng)建一個連續(xù)分布，繼承rv_continuous類創(chuàng)建連續(xù)分布是非常簡單的。>>>>>>fromscipyimportstats>>>classdeterministic_gen(stats.rv_continuous):.def_cdf(self,x):.returnnp.where(x<0,

25、0.,1.).def_stats(self):.return0.,0.,0.,0.>>>>>>deterministic=deterministic_gen(name="deterministic")>>>deterministic.cdf(np.arange(-3,3,0.5)array(0.,0.,0.,0.,0.,0.,1.,1.,1.,1.,1.,1.)令人高興的是，pdf也能被自動計算出來：>>>>>>deterministic.pdf(np.arange(-3,3,0.5)

26、array( 0.00000000e+00,0.00000000e+00,0.00000000e+00,0.00000000e+00,0.00000000e+00,0.00000000e+00,5.83333333e+04,4.16333634e-12,4.16333634e-12,4.16333634e-12,4.16333634e-12,4.16333634e-12)注意這種用法的性能問題,參見“性能問題與注意事項”一節(jié)。這種缺乏信息的通用計算可能非常慢。而且再看看下面這個準確性的例子:>>>>>>egrateimportqua

27、d>>>quad(deterministic.pdf,-1e-1,1e-1)(4.163336342344337e-13,0.0)但這并不是對pdf積分的正確的結(jié)果，實際上結(jié)果應為1.讓我們將積分變得更小一些。>>>>>>quad(deterministic.pdf,-1e-3,1e-3)#warningremoved(1.0003,0.18182458)這樣看上去好多了，然而，問題本身來源于pdf不是來自包給定的類的定義。繼承rv_discrete類在之后我們使用stats.rv_discrete產(chǎn)生一個離散分布，其有一個整數(shù)區(qū)間截斷概率

28、。通用信息通用信息可以從rv_discrete的docstring中得到>>>>>>fromscipy.statsimportrv_discrete>>>help(rv_discrete)我們學到這一點：“你可以構建任意一個像P(X=xk尸pk一樣形式的離散rv,通過彳專遞(xk,pk)元組序列給rv_discrete初始化方法(通過value=keyword方式)，但其不能有0概率值?！苯酉聛?，還有一些進一步的要求：關鍵字名字必須知道。Xk點必須是整數(shù)小數(shù)的有效位數(shù)應當被給出。事實上，如果最后兩個要求沒有被滿足，一個異常將被拋出或者導致一

29、個錯誤的數(shù)值。一個例子讓我們開始辦，首先>>>> >>npoints=20#numberofintegersupportpointsofthedistributionminus1> >>npointsh=npoints/2> >>npointsf=float(npoints)> >>nbound=4#boundsforthetruncatednormal> >>normbound=(1+1/npointsf)*nbound#actualboundsoftruncatednormal>

30、; >>grid=np.arange(-npointsh,npointsh+2,1)#integergrid> >>gridlimitsnorm=(grid-0.5)/npointsh*nbound#binlimitsforthetruncnorm> >>gridlimits=grid-0.5#usedlaterintheanalysis> >>grid=grid:-1> >>probs=np.diff(stats.truncnorm.cdf(gridlimitsnorm,-normbound,normboun

31、d)> >>gridint=grid最后我們可以繼承rv_discrete類>>>> >>normdiscrete=stats.rv_discrete(values=(gridint,.np.round(probs,decimals=7),name='normdiscrete')現(xiàn)在我們已經(jīng)定義了這個分布，我們可以調(diào)用其所有常規(guī)的離散分布方法。>>>> >>print'mean=%6.4f,variance=%6.4f,skew=%6.4f,kurtosis=%6.4f'%

32、.normdiscrete.stats(moments='mvsk')mean=-0.0000,variance=6.3302,skew=0.0000,kurtosis=-0.0076>>>> >>nd_std=np.sqrt(normdiscrete.stats(moments='v')測試上面的結(jié)果讓我們產(chǎn)生一個隨機樣本并且比較連續(xù)概率的情況。>>>> >>n_sample=500>>>np.random.seed()#fixtheseedforreplicabilit

33、y>>>rvs=normdiscrete.rvs(size=n_sample)>>>rvsnd=rvs>>>f,l=np.histogram(rvs,bins=gridlimits)>>>sfreq=np.vstack(gridint,f,probs*n_sample).T>>>printsfreq-1.00000000e+010.00000000e+00-9.00000000e+000.00000000e+00-8.00000000e+000.00000000e+00-7.00000000e+002.0

34、0000000e+00-6.00000000e+001.00000000e+00-5.00000000e+009.00000000e+00-4.00000000e+002.60000000e+01-3.00000000e+003.70000000e+01-2.00000000e+005.10000000e+01-1.00000000e+007.10000000e+010.00000000e+007.40000000e+011.00000000e+008.90000000e+012.00000000e+005.50000000e+013.00000000e+005.00000000e+014.0

35、0000000e+001.70000000e+015.00000000e+001.10000000e+012.95019349e-021.32294142e-015.06497902e-011.65568919e+004.62125309e+001.10137298e+012.24137683e+013.89503370e+015.78004747e+017.32455414e+017.92618251e+017.32455414e+015.78004747e+013.89503370e+012.24137683e+011.10137298e+016.00000000e+004.0000000

36、0e+007.00000000e+003.00000000e+008.00000000e+000.00000000e+009.00000000e+000.00000000e+001.00000000e+010.00000000e+004.62125309e+001.65568919e+005.06497902e-011.32294142e-012.95019349e-02(Sourcecode)-IU9-7-6.540123457a910(Sourcecode)接下來，我們可以測試，是否我們的樣本取自于一個normdiscrete分布。這也是在驗證是否隨機數(shù)是以正確的方式產(chǎn)生的?？ǚ綔y試要求起

37、碼在每個子區(qū)間(bin)里具有最小數(shù)目的觀測值。我們組合末端子區(qū)間進大子區(qū)間所以它們現(xiàn)在包含了足夠數(shù)量的觀測值。>>>>>>f2=np.hstack(f:5.sum(),f5:-5,f-5:.sum()>>>p2=np.hstack(probs:5.sum(),probs5:-5,probs-5:.sum()>>>ch2,pval=stats.chisquare(f2,p2*n_sample)>>>>>>print'chisquarefornormdiscrete:chi2=%6

38、.3fpvalue=%6.4f'%(ch2,pval)chisquarefornormdiscrete:chi2=12.466pvalue=0.4090P值在這個情況下是不顯著地，所以我們可以斷言我們的隨機樣本的確是由此分布產(chǎn)生的。樣本分析首先，我們創(chuàng)建一些隨機變量。我們設置一個種子所以每次我們都可以得到相同的結(jié)果以便觀察。作為一個例子，我們從t分布中抽一個樣本。>>>>>>np.random.seed(282629734)>>>x=stats.t.rvs(10,size=1000)這里，我們設置了t分布的形態(tài)參數(shù)，在這里就是自由度

39、，設為10.使用size=1000表示我們的樣本由1000個抽樣是獨立的(偽)。當我們不指派10c和scale時，它們具有默認值0和1.描述統(tǒng)計X是一個numpy數(shù)組。我們可以直接調(diào)用它所有的方法。>>>> >>printx.max(),x.min()#equivalenttonp.max(x),np.min(x)5.26327732981-3.78975572422> >>printx.mean(),x.var()#equivalenttonp.mean(x),np.var(x)0.851.28899386208如何比較分布本身和它的樣本

40、的指標？>>>> >>m,v,s,k=stats.t.stats(10,moments='mvsk')>>>n,(smin,smax),sm,sv,ss,sk=stats.describe(x)>>>>>>print'distribution:',distribution:>>>sstr='mean=%6.4f,variance=%6.4f,skew=%6.4f,kurtosis=%6.4f'>>>printsstr%(m

41、,v,s,k)mean=0.0000,variance=1.2500,skew=0.0000,kurtosis=1.0000>>>print'sample:',sample:>>>printsstr%(sm,sv,ss,sk)mean=0.0141,variance=1.2903,skew=0.2165,kurtosis=1.0556注意：stats.describe用的是無偏的方差估計量，而np.var卻用的是有偏的估計量。T檢驗和KS檢驗我們可以使用t檢驗是否樣本與給定均值(這里是理論均值)存在統(tǒng)計顯著差異。>>>>

42、; >>print't-statistic=%6.3fpvalue=%6.4f'%stats.ttest_1samp(x,m)t-statistic=0.391pvalue=0.6955P值為0.7,這代表第一類錯誤的概率，在例子中，為10%o我們不能拒絕“該樣本均值為0”這個假設，0是標準t分布的理論均值。>>>> >>tt=(sm-m)/np.sqrt(sv/float(n)#t-statisticformean> >>pval=stats.t.sf(np.abs(tt),n-1)*2#two-sidedpv

43、alue=Prob(abs(t)>tt)> >>print't-statistic=%6.3fpvalue=%6.4f'%(tt,pval)t-statistic=0.391pvalue=0.6955這里Kolmogorov-Smirnov檢驗(KS檢驗)被被用來檢驗樣本是否來自一個標準的t分布。>>>>>>print'KS-statisticD=%6.3fpvalue=%6.4f'%stats.kstest(x,'t',(10,)KS-statisticD=0.016pvalue=0.

44、9606又一次得到了很高的P值。所以我們不能拒絕樣本是來自t分布的假設。在實際應用中，我們不能知道潛在的分布到底是什么。如果我們使用KS檢驗我們的樣本對照正態(tài)分布，那么我們將也不能拒絕我們的樣本是來自正態(tài)分布的，在這種情況下P值為0.4左右。>>>>>>print'KS-statisticD=%6.3fpvalue=%6.4f'%stats.kstest(x,'norm')KS-statisticD=0.028pvalue=0.3949無論如何，標準正態(tài)分布有1的方差，當我們的樣本有1.29時。如果我們標準化我們的樣本并且測試

45、它比照正態(tài)分布，那么P值將又一次很高我們將還是不能拒絕假設是來自正態(tài)分布的。>>>>>>d,pval=stats.kstest(x-x.mean()/x.std(),'norm')>>>print'KS-statisticD=%6.3fpvalue=%6.4f'%(d,pval)KS-statisticD=0.032pvalue=0.2402注釋：KS檢驗假設我們比照的分布就是以給定的參數(shù)確定的，但我們在最后估計了均值和方差，這個假設就被違反了，故而這個測試統(tǒng)計量的P值是含偏的，這個用法是錯誤的。分布尾部最后

46、，我們可以檢查分布的右尾，我們可以使用分位點函數(shù)ppf,其為cdf函數(shù)的逆，來獲得臨界值,或者更直接的，我們可以使用殘存函數(shù)的逆來辦。>>>>>> crit01, crit05,crit10 = stats .t.ppf(1-0.01 , 1-0.05 , 1-0.10, 10)>>> print 'criticalvalues fromppfat1%, 5% and10% %8.4f %8.4f %8.4f'% (crit01, crit05,crit10)critical values fromppf at 1%,5%a

47、nd10%2.76381.81251.3722>>> print 'criticalvalues fromisfat 1%,5% and10% %8.4f %8.4f %8.4f'% tuple (stats.t.isf(0.01,0.05,0.10,10)critical values from isf at 1%, 5%and10%2.76381.81251.3722>>>>>>freq01 = np .sum(x > crit01)float100>>>freq05 = np .sum(x &g

48、t; crit05)float100>>>freq10 = np .sum(x > crit10)float100>>>print 'sample %-frequencyat 1%,5% and10% tail%8.4f %8.4f %8.4f'% (freq01,freq05,freq10)sample %-frequency at 1%, 5% and 10% tail 1.40005.800010.5000在這三種情況中，我們的樣本有有一個更重的尾部，即實際在理論分界值右邊的概率要高于理論值。我們可以通過使用更大的樣本來獲得更好的

49、擬合。在以下情況經(jīng)驗頻率已經(jīng)很接近理論概率了，但即使我們重復這個過程若干次，波動依然會保持在這個程度。>>>>>>freq05l=np.sum(stats.t.rvs(10,size=10000)>crit05)/10000.0*100>>>print'largersample%-frequencyat5%tail%8.4f'%freq05llargersample%-frequencyat5%tail4.8000我們也可以比較它與正態(tài)分布的尾部，其有一個輕的多的尾部：>>>> >>

50、print'tailprob.ofnormalat1%,5%and10%8.4f%8.4f%8.4f'%.tuple(stats.norm.sf(crit01,crit05,crit10)*100)tailprob.ofnormalat1%,5%and10%0.28573.49578.5003卡方檢驗可以被用來測試，是否一個有限的分類觀測值頻率與假定的理論概率分布具有顯著差異。>>>> >>quantiles=0.0,0.01,0.05,0.1,1-0.10,1-0.05,1-0.01,1.0> >>crit=stats.t

51、.ppf(quantiles,10)> >>printcrit-Inf-2.76376946-1.81246112-1.372183641.372183641.812461122.76376946Inf> >>n_sample=x.size> >>freqcount=np.histogram(x,bins=crit)0> >>tprob=np.diff(quantiles)> >>nprob=np.diff(stats.norm.cdf(crit)> >>tch,tpval=stats.

52、chisquare(freqcount,tprob*n_sample)> >>nch,npval=stats.chisquare(freqcount,nprob*n_sample)> >>print'chisquarefort:chi2=%6.3fpvalue=%6.4f'%(tch,tpval)chisquare for t:chi22.300 pvalue0.8901> >>print'chisquarefornormal:chi2=%6.3fpvalue=%6.4f'%(nch,npval)chisqu

53、arefornormal:chi2=64.605pvalue=0.0000我們看到當t分布檢驗沒被拒絕時標準正態(tài)分布卻被完全拒絕。在我們的樣本區(qū)分出這兩個分布后，我們可以先進行擬合確定scale與location再檢查擬合后的分布的差異性。我們可以先進行擬合，再用擬合分布而不是默認(起碼location和scale是默認的)分布去進行檢驗。>>>> >>tdof,tloc,tscale=stats.t.fit(x)> >>nloc,nscale=stats.norm.fit(x)> >>tprob=np.diff(stat

54、s.t.cdf(crit,tdof,loc=tloc,scale=tscale)> >>nprob=np.diff(stats.norm.cdf(crit,loc=nloc,scale=nscale)> >>tch,tpval=stats.chisquare(freqcount,tprob*n_sample)> >>nch,npval=stats.chisquare(freqcount,nprob*n_sample)> >>print'chisquarefort:chi2=%6.3fpvalue=%6.4f'

55、;%(tch,tpval)chisquarefort:chi2=1.577pvalue=0.9542> >>print'chisquarefornormal:chi2=%6.3fpvalue=%6.4f'%(nch,npval)chisquarefornormal:chi2=11.084pvalue=0.0858在經(jīng)過參數(shù)調(diào)整之后，我們?nèi)匀豢梢砸?%水平拒絕正態(tài)分布假設。然而卻以95%的p值顯然的不能拒絕t分布。正態(tài)分布的特殊檢驗自從正態(tài)分布變?yōu)榻y(tǒng)計學中最常見的分布，就出現(xiàn)了大量的方法用來檢驗一個樣本是否可以被看成是來自正態(tài)分布的。首先我們檢驗分布的峰度和偏度

56、是否顯著地與正態(tài)分布的對應值相差異。>>>> >>print'normalskewtestteststat=%6.3fpvalue=%6.4f'%stats.skewtest(x)normalskewtestteststat=2.785pvalue=0.0054> >>print'normalkurtosistestteststat=%6.3fpvalue=%6.4f'%stats.kurtosistest(x)normalkurtosistestteststat=4.757pvalue=0.0000將這兩個

57、檢驗組合起來的正態(tài)性檢驗>>>> >>print'normaltestteststat=%6.3fpvalue=%6.4f'%stats.normaltest(x)normaltestteststat=30.379pvalue=0.0000在所有三個測試中，P值是非常低的，所以我們可以拒絕我們的樣本的峰度與偏度與正態(tài)分布相同的假設。當我們的樣本標準化之后，我們依舊得到相同的結(jié)果。>>>> >>print'normaltestteststat=%6.3fpvalue=%6.4f'%.stats

58、.normaltest(x-x.mean()/x.std()normaltestteststat=30.379pvalue=0.0000因為正態(tài)性被很強的拒絕了，所以我們可以檢查這種檢驗方式是否可以有效地作用到其他情況中。>>>> >>print'normaltestteststat=%6.3fpvalue=%6.4f'%stats.normaltest(stats.t.rvs(10,size=100)normaltestteststat=4.698pvalue=0.0955> >>print'normaltestt

59、eststat=%6.3fpvalue=%6.4f'%stats.normaltest(stats.norm.rvs(size=1000)normaltestteststat=0.613pvalue=0.7361我們檢驗了小樣本的t分布樣本的觀測值以及一個大樣本的正態(tài)分布觀測值，在這兩種情況中我們都不能拒絕其來自正態(tài)分布的空假設。得到這樣的結(jié)果是因為前者是因為無法區(qū)分小樣本下的t分布,后者是因為它本來就來自正態(tài)分布。比較兩個樣本接下來，我們有兩個分布，其可以判定為相同或者來自不同的分布，以及我們希望測試是否這些樣本有相同的統(tǒng)計特征。均值c:iknowdocsharedatacurwor

60、kncetutorialstats.html以相同的均值產(chǎn)生的樣本進行檢驗：>>>> >>rvs1=stats.norm.rvs(loc=5,scale=10,size=500)> >>rvs2=stats.norm.rvs(loc=5,scale=10,size=500)> >>stats.ttest_ind(rvs1,rvs2)(-0.548983,0.5831943748663857)以不同的均值產(chǎn)生的樣本進行檢驗：>>>> >>rvs3=stats.norm.rvs(loc=8,

61、scale=10,size=500)> >>stats.ttest_ind(rvs1,rvs3)(-4.5334142901750321,6.5895e-006)對于兩個不同的樣本進行的KS檢驗在這個例子中我們使用兩個同分布的樣本進行檢驗.設因為p值很高，毫不奇怪我們不能拒絕原假。>>>> >>stats.ks_2samp(rvs1,rvs2)(0.9999995,0.99541195173064878)在第二個例子中，由于均值不同，所以我們可以拒絕空假設，由P值小于1%。>>>> >>stats.ks_

62、2samp(rvs1,rvs3)(0.199999,0.00273141)核密度估計一個常見的統(tǒng)計學問題是從一個樣本中估計隨機變量的概率密度分布函數(shù)(PDF)這個問題被稱為密度估計，對此最著名的工具是直方圖。直方圖是一個很好的可視化工具(主要是因為每個人都理解它)但是對于對于數(shù)據(jù)特征的利用卻并不是非常有效率。核密度估計(KDE對于這個問題)是一個更有效的工具。這個gaussian_kde估計方法可以被用來估at單元或多元數(shù)據(jù)的PDF。它在數(shù)據(jù)呈單峰的時候工作的最好，但也可以在多峰情況下工作。單元估計c:iknowdocsharedatacurworkyreferencetutorialstats.html我們以一個最小數(shù)據(jù)集來觀察gaussian_kde是如何工作的，以及帶寬(bandwidth)的不同選擇方式。PDF對應的數(shù)據(jù)被以藍線的形式顯示在圖像的底端(被稱為毯圖(rugplot)>>>>>>fromscipyimportstats> >>importmatplotlib.pyplotasplt>>>> >>x1=np.array(-7,-5,1,4,5,dtype=np.float)> >>kde1=st