Python統(tǒng)計學(xué)包scipy.stats手冊

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

介紹在這個教程我們討論一部分scipy.stats模塊的特性。這里我們的意圖是提供應(yīng)使用者一個關(guān)于這個包的實用性知識。我們推薦referencemanual來介紹更多的細節(jié)。注意：這個文檔還在發(fā)展中。隨機變量有一些通用的概率分布類被封裝在continuousrandomvariables以及discreterandomvariables中。有80多個連續(xù)性隨機變量(RVs)以及10余個離散隨機變量已經(jīng)用這些類建立。同樣，新的程序和分布可以被用戶新建（假如你構(gòu)造了一個，請?zhí)峁┧o我們幫助發(fā)展這個包）。所有記錄函數(shù)被放在子包scipy.stats中，且有這些函數(shù)的一個幾乎完整的列表可以使用info(stats)獲得。這個列表里的隨機變量也可以從stats子包的docstring中獲得介紹。在接下來的討論中，我們著重于連續(xù)性隨機變量(RVs)。幾乎所有離散變量也符合下面的討論，但是我們也要指出一些區(qū)別在“離散分布的特殊之處”中。獲得幫助所有分布可以使用help函數(shù)得到解釋。為獲得這些信息只需要使用像這樣的簡樸調(diào)用：>>>>>>fromscipyimportstats>>>fromscipy.statsimportnorm>>>printnorm.__doc__作為例子，我們用這種方式找分布的上下界>>>>>>print'boundsofdistributionlower:%s,upper:%s'%(norm.a,norm.b)boundsofdistributionlower:-inf,upper:inf我們可以通過調(diào)用dir(norm)來獲得關(guān)于這個（正態(tài)）分布的所有方法和屬性。應(yīng)當(dāng)看到，一些方法是私有方法盡管其并沒有以名稱表達出來（比如它們前面沒有以下劃線開頭），比如veccdf就只用于內(nèi)部計算（試圖使用那些方法將引發(fā)警告，由于它們也許會在后續(xù)開發(fā)中被移除）為了獲得真正的重要方法，我們列舉凍結(jié)分布的方法（我們將在下文解釋何謂“凍結(jié)分布”）>>>>>>rv=norm()>>>dir(rv)#reformatted['__class__','__delattr__','__dict__','__doc__','__getattribute__','__hash__','__init__','__module__','__new__','__reduce__','__reduce_ex__','__repr__','__setattr__','__str__','__weakref__','args','cdf','dist','entropy','isf','kwds','moment','pdf','pmf','ppf','rvs','sf','stats']最后，我們能通過內(nèi)省獲得所有的可用分布的信息。>>>>>>importwarnings>>>warnings.simplefilter('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_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)，（費舍爾）峰度。moment:分布的非中心矩。讓我們使用一個標準的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])相應(yīng)的，像pdf,cdf之類的簡樸方法可以用np.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)array([-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")(array(3.0),array(16.0))通常經(jīng)標準化的分布的隨機變量X可以通過變換(X-loc)/scale獲得。它們的默認值是loc=0以及scale=1.聰明的使用loc與scale可以幫助以靈活的方式調(diào)整標準分布達成所想要的效果。為了進一步說明縮放的效果，下面給出盼望為1/λ指數(shù)分布的cdf。F(x)=1?exp(?λx)通過像上面那樣使用scale，可以看到得到想要的盼望值。>>>>>>fromscipy.statsimportexpon>>>expon.mean(scale=3.)3.0均勻分布也是令人感愛好的：>>>>>>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.704在這里，為解釋最后一段的輸出：norm.rvs(5)產(chǎn)生了一個正態(tài)分布變量，其盼望，即loc=5.我傾向于明確的使用loc,scale作為關(guān)鍵字而非像上面那樣運用參數(shù)的順序。由于這看上去有點令人困惑。我們在我們解釋“凍結(jié)RV”的主題之前澄清這一點。形態(tài)參數(shù)雖然一個一般的連續(xù)隨機變量可以通過賦予loc和scale參數(shù)擬定，但一些分布還需要額外的形態(tài)參數(shù)。作為例子，看到這個伽馬分布，這是它的密度函數(shù)γ(x,a)=λ(λx)a?1Γ(a)e?λx,規(guī)定一個形態(tài)參數(shù)a。注意到λ的設(shè)立可以通過設(shè)立scale關(guān)鍵字為1/λ進行。讓我們檢查伽馬分布的形態(tài)參數(shù)的名字的數(shù)量。（我們知道從上面知道其應(yīng)當(dāng)為1）>>>>>>fromscipy.statsimportgamma>>>gamma.numargs1>>>gamma.shapes'a'現(xiàn)在我們設(shè)立形態(tài)變量的值為1以變成指數(shù)分布。所以我們可以容易的比較是否得到了我們所盼望的結(jié)果。>>>>>>gamma(1,scale=2.).stats(moments="mv")(array(2.0),array(4.0))注意我們也可以以關(guān)鍵字的方式指定形態(tài)參數(shù)：>>>>>>gamma(a=1,scale=2.).stats(moments="mv")(array(2.0),array(4.0))凍結(jié)分布不斷地傳遞loc與scale關(guān)鍵字最終會讓人厭煩。而凍結(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)這正是我們應(yīng)當(dāng)?shù)玫降?。廣播像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)用了兩次isf產(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)，則進行相應(yīng)匹配，我們可以分別得到10%，5%，1%尾的臨界值對于10，11,12的自由度。>>>>>>stats.t.isf([0.1,0.05,0.01],[10,11,12])array([1.37218364,1.79588482,2.68099799])離散分布的特殊之處離散分布的方法的大多數(shù)與連續(xù)分布很類似。當(dāng)然像pdf被更換為密度函數(shù)pmf，沒有估計方法，像fit就不能用了。而scale不是一個合法的關(guān)鍵字參數(shù)。Location參數(shù)，關(guān)鍵字loc則仍然可以使用用于位移。cdf的計算規(guī)定一些額外的關(guān)注。在連續(xù)分布的情況下，累積分布函數(shù)在大多數(shù)標準情況下是嚴格遞增的，所以有唯一的逆。而cdf在離散分布，無論如何，是階躍函數(shù)，所以cdf的逆，分位點函數(shù)，規(guī)定一個不同的定義：ppf(q)=min{x: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=hypergeom.cdf(x,M,n,N)>>>prbarray([0.4066,0.05211,0.9301,0.7386])>>>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.])分布擬合非凍結(jié)分布的參數(shù)估計的重要方法：fit：分布參數(shù)的極大似然估計，涉及l(fā)ocation與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的rvs特殊函數(shù)給出。一般來說使用精確計算會比較快。另一方面，通用方法被用于當(dāng)分布沒有被指派明確的計算方法時使用。為了定義一個分布，只有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ū)間上的值被測試過了，然而在一些奇葩的邊界，仍然也許有錯誤的值存在。fit的極大似然估計以默認值作為初始參數(shù)將不會工作的很好，用戶必須指派合適的初始參數(shù)。并且，對于一些分布使用極大似然估計自身就不是一個好的選擇。構(gòu)造具體的分布下一個例子展示了如何建立你自己的分布。更多的例子見分布用法以及記錄檢查創(chuàng)建一個連續(xù)分布，繼承rv_continuous類創(chuàng)建連續(xù)分布是非常簡樸的。>>>>>>fromscipyimportstats>>>classdeterministic_gen(stats.rv_continuous):...def_cdf(self,x):...returnnp.where(x<0,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))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é)。這種缺少信息的通用計算也許非常慢。并且再看看下面這個準確性的例子：>>>>>>egrateimportquad>>>quad(deterministic.pdf,-1e-1,1e-1)(4.4337e-13,0.0)但這并不是對pdf積分的對的的結(jié)果，事實上結(jié)果應(yīng)為1.讓我們將積分變得更小一些。>>>>>>quad(deterministic.pdf,-1e-3,1e-3)#warningremoved(1.173,0.8182458)這樣看上去好多了，然而，問題自身來源于pdf不是來自包給定的類的定義。繼承rv_discrete類在之后我們使用stats.rv_discrete產(chǎn)生一個離散分布，其有一個整數(shù)區(qū)間截斷概率。通用信息通用信息可以從rv_discrete的docstring中得到。>>>>>>fromscipy.statsimportrv_discrete>>>help(rv_discrete)我們學(xué)到這一點：“你可以構(gòu)建任意一個像P(X=xk)=pk同樣形式的離散rv，通過傳遞(xk,pk)元組序列給rv_discrete初始化方法(通過value=keyword方式)，但其不能有0概率值?！苯酉聛恚杏幸恍┻M一步的規(guī)定：關(guān)鍵字名字必須知道。Xk點必須是整數(shù)小數(shù)的有效位數(shù)應(yīng)當(dāng)被給出。事實上，假如最后兩個規(guī)定沒有被滿足，一個異常將被拋出或者導(dǎo)致一個錯誤的數(shù)值。一個例子讓我們開始辦，一方面>>>>>>npoints=20#numberofintegersupportpointsofthedistributionminus1>>>npointsh=npoints/2>>>npointsf=float(npoints)>>>nbound=4#boundsforthetruncatednormal>>>normbound=(1+1/npointsf)*nbound#actualboundsoftruncatednormal>>>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,normbound))>>>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'%\...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(87655678)#fixtheseedforreplicability>>>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+002.95019349e-02][-9.00000000e+000.00000000e+001.32294142e-01][-8.00000000e+000.00000000e+005.06497902e-01][-7.00000000e+002.00000000e+001.65568919e+00][-6.00000000e+001.00000000e+004.62125309e+00][-5.00000000e+009.00000000e+001.10137298e+01][-4.00000000e+002.60000000e+012.24137683e+01][-3.00000000e+003.70000000e+013.89503370e+01][-2.00000000e+005.10000000e+015.78004747e+01][-1.00000000e+007.10000000e+017.32455414e+01][0.00000000e+007.40000000e+017.92618251e+01][1.00000000e+008.90000000e+017.32455414e+01][2.00000000e+005.50000000e+015.78004747e+01][3.00000000e+005.00000000e+013.89503370e+01][4.00000000e+001.70000000e+012.24137683e+01][5.00000000e+001.10000000e+011.10137298e+01][6.00000000e+004.00000000e+004.62125309e+00][7.00000000e+003.00000000e+001.65568919e+00][8.00000000e+000.00000000e+005.06497902e-01][9.00000000e+000.00000000e+001.32294142e-01][1.00000000e+010.00000000e+002.95019349e-02]](Sourcecode)(Sourcecode)接下來，我們可以測試，是否我們的樣本取自于一個normdiscrete分布。這也是在驗證是否隨機數(shù)是以對的的方式產(chǎn)生的。卡方測試規(guī)定起碼在每個子區(qū)間（bin）里具有最小數(shù)目的觀測值。我們組合末端子區(qū)間進大子區(qū)間所以它們現(xiàn)在包含了足夠數(shù)量的觀測值。>>>>>>f2=np.hstack([f[:5].sum(),f[5:-5],f[-5:].sum()])>>>p2=np.hstack([probs[:5].sum(),probs[5:-5],probs[-5:].sum()])>>>ch2,pval=stats.chisquare(f2,p2*n_sample)>>>>>>print'chisquarefornormdiscrete:chi2=%6.3fpvalue=%6.4f'%(ch2,pval)chisquarefornormdiscrete:chi2=12.466pvalue=0.4090P值在這個情況下是不顯著地，所以我們可以斷言我們的隨機樣本的確是由此分布產(chǎn)生的。樣本分析一方面，我們創(chuàng)建一些隨機變量。我們設(shè)立一個種子所以每次我們都可以得到相同的結(jié)果以便觀測。作為一個例子，我們從t分布中抽一個樣本。>>>>>>np.random.seed()>>>x=stats.t.rvs(10,size=1000)這里，我們設(shè)立了t分布的形態(tài)參數(shù)，在這里就是自由度，設(shè)為10.使用size=1000表達我們的樣本由1000個抽樣是獨立的（偽）。當(dāng)我們不指派loc和scale時，它們具有默認值0和1.描述記錄X是一個numpy數(shù)組。我們可以直接調(diào)用它所有的方法。>>>>>>printx.max(),x.min()#equivalenttonp.max(x),np.min(x)5.-3.>>>printx.mean(),x.var()#equivalenttonp.mean(x),np.var(x)0.051.如何比較分布自身和它的樣本的指標？>>>>>>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,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檢查是否樣本與給定均值（這里是理論均值）存在記錄顯著差異。>>>>>>print't-statistic=%6.3fpvalue=%6.4f'%stats.ttest_1samp(x,m)t-statistic=0.391pvalue=0.6955P值為0.7，這代表第一類錯誤的概率，在例子中，為10%。我們不能拒絕“該樣本均值為0”這個假設(shè)，0是標準t分布的理論均值。>>>>>>tt=(sm-m)/np.sqrt(sv/float(n))#t-statisticformean>>>pval=stats.t.sf(np.abs(tt),n-1)*2#two-sidedpvalue=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.9606又一次得到了很高的P值。所以我們不能拒絕樣本是來自t分布的假設(shè)。在實際應(yīng)用中，我們不能知道潛在的分布到底是什么。假如我們使用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的方差，當(dāng)我們的樣本有1.29時。假如我們標準化我們的樣本并且測試它比照正態(tài)分布，那么P值將又一次很高我們將還是不能拒絕假設(shè)是來自正態(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è)我們比照的分布就是以給定的參數(shù)擬定的，但我們在最后估計了均值和方差，這個假設(shè)就被違反了，故而這個測試記錄量的P值是含偏的，這個用法是錯誤的。分布尾部最后，我們可以檢查分布的右尾，我們可以使用分位點函數(shù)ppf，其為cdf函數(shù)的逆，來獲得臨界值，或者更直接的，我們可以使用殘存函數(shù)的逆來辦。>>>>>>crit01,crit05,crit10=stats.t.ppf([1-0.01,1-0.05,1-0.10],10)>>>print'criticalvaluesfromppfat1%%,5%%and10%%%8.4f%8.4f%8.4f'%(crit01,crit05,crit10)criticalvaluesfromppfat1%,5%and10%2.76381.81251.3722>>>print'criticalvaluesfromisfat1%%,5%%and10%%%8.4f%8.4f%8.4f'%tuple(stats.t.isf([0.01,0.05,0.10],10))criticalvaluesfromisfat1%,5%and10%2.76381.81251.3722>>>>>>freq01=np.sum(x>crit01)/float(n)*100>>>freq05=np.sum(x>crit05)/float(n)*100>>>freq10=np.sum(x>crit10)/float(n)*100>>>print'sample%%-frequencyat1%%,5%%and10%%tail%8.4f%8.4f%8.4f'%(freq01,freq05,freq10)sample%-frequencyat1%,5%and10%tail1.40005.800010.5000在這三種情況中，我們的樣本有有一個更重的尾部，即實際在理論分界值右邊的概率要高于理論值。我們可以通過使用更大的樣本來獲得更好的擬合。在以下情況經(jīng)驗頻率已經(jīng)很接近理論概率了，但即使我們反復(fù)這個過程若干次，波動仍然會保持在這個限度。>>>>>>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)分布的尾部，其有一個輕的多的尾部：>>>>>>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.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.chisquare(freqcount,tprob*n_sample)>>>nch,npval=stats.chisquare(freqcount,nprob*n_sample)>>>print'chisquarefort:chi2=%6.3fpvalue=%6.4f'%(tch,tpval)chisquarefort:chi2=2.300pvalue=0.8901>>>print'chisquarefornormal:chi2=%6.3fpvalue=%6.4f'%(nch,npval)chisquarefornormal:chi2=64.605pvalue=0.0000我們看到當(dāng)t分布檢查沒被拒絕時標準正態(tài)分布卻被完全拒絕。在我們的樣本區(qū)分出這兩個分布后，我們可以先進行擬合擬定scale與location再檢查擬合后的分布的差異性。我們可以先進行擬合，再用擬合分布而不是默認（起碼location和scale是默認的）分布去進行檢查。>>>>>>tdof,tloc,tscale=stats.t.fit(x)>>>nloc,nscale=stats.norm.fit(x)>>>tprob=np.diff(stats.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'%(tch,tpval)chisquarefort:chi2=1.577pvalue=0.9542>>>print'chisquarefornormal:chi2=%6.3fpvalue=%6.4f'%(nch,npval)chisquarefornormal:chi2=11.084pvalue=0.0858在通過參數(shù)調(diào)整之后，我們?nèi)匀豢梢砸?%水平拒絕正態(tài)分布假設(shè)。然而卻以95%的p值顯然的不能拒絕t分布。正態(tài)分布的特殊檢查自從正態(tài)分布變?yōu)橛涗泴W(xué)中最常見的分布，就出現(xiàn)了大量的方法用來檢查一個樣本是否可以被當(dāng)作是來自正態(tài)分布的。一方面我們檢查分布的峰度和偏度是否顯著地與正態(tài)分布的相應(yīng)值相差異。>>>>>>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將這兩個檢查組合起來的正態(tài)性檢查>>>>>>print'normaltestteststat=%6.3fpvalue=%6.4f'%stats.normaltest(x)normaltestteststat=30.379pvalue=0.0000在所有三個測試中，P值是非常低的，所以我們可以拒絕我們的樣本的峰度與偏度與正態(tài)分布相同的假設(shè)。當(dāng)我們的樣本標準化之后，我們依舊得到相同的結(jié)果。>>>>>>print'normaltestteststat=%6.3fpvalue=%6.4f'%\...stats.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'normaltestteststat=%6.3fpvalue=%6.4f'%stats.normaltest(stats.norm.rvs(size=1000))normaltestteststat=0.613pvalue=0.7361我們檢查了小樣本的t分布樣本的觀測值以及一個大樣本的正態(tài)分布觀測值，在這兩種情況中我們都不能拒絕其來自正態(tài)分布的空假設(shè)。得到這樣的結(jié)果是由于前者是由于無法區(qū)分小樣本下的t分布，后者是由于它本來就來自正態(tài)分布。比較兩個樣本接下來，我們有兩個分布，其可以鑒定為相同或者來自不同的分布，以及我們希望測試是否這些樣本有相同的記錄特性。均值以相同的均值產(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.88583,0.58357)以不同的均值產(chǎn)生的樣本進行檢查：>>>>>>rvs3=stats.norm.rvs(loc=8,scale=10,size=500)>>>stats.ttest_ind(rvs1,rvs3)(-4.53341,6.895e-006)對于兩個不同的樣本進行的KS檢查在這個例子中我們使用兩個同分布的樣本進行檢查.設(shè)由于P值很高,毫不奇怪我們不能拒絕原假。>>>>>>stats.ks_2samp(rvs1,rvs2)(0.999995,0.995418)在第二個例子中，由于均值不同，所以我們可以拒絕空假設(shè)，由P值小于1%。>>>>>>stats.ks_2samp(rvs1,rvs3)(0.199999,0.00273141)核密度估計一個常見的記錄學(xué)問題是從一個樣本中估計隨機變量的概率密度分布函數(shù)（PDF）這個問題被稱為密度估計，對此最著名的工具是直方圖。直方圖是一個很好的可視化工具（重要是由于每個人都理解它）。但是對于對于數(shù)據(jù)特性的運用卻并不是非常有效率。核密度估計(KDE對于這個問題)是一個更有效的工具。這個gaussian_kde估計方法可以被用來估計單元或多元數(shù)據(jù)的PDF。它在數(shù)據(jù)呈單峰的時候工作的最佳，但也可以在多峰情況下工作。單元估計我們以一個最小數(shù)據(jù)集來觀測gaussian_kde是如何工作的，以及帶寬（bandwidth）的不同選擇方式。PDF相應(yīng)的數(shù)據(jù)被以藍線的形式顯示在圖像的底端（被稱為毯圖（rugplot））>>>>>>fromscipyimportstats>>>importmatplotlib.pyplotasplt>>>>>>x1=np.array([-7,-5,1,4,5],dtype=np.float)>>>kde1=stats.gaussian_kde(x1)>>>kde2=stats.gaussian_kde(x1,bw_method='silverman')>>>>>>fig=plt.figure()>>>ax=fig.add_subplot(111)>>>>>>ax.plot(x1,np.zeros(x1.shape),'b+',ms=20)#rugplot>>>x_eval=np.linspace(-10,10,num=200)>>>ax.plot(x_eval,kde1(x_eval),'k-',label="Scott'sRule")>>>ax.plot(x_eval,kde1(x_eval),'r-',label="Silverman'sRule")>>>>>>plt.show()(Sourcecode)我們看到在Scott規(guī)則以及Silverman規(guī)則下的結(jié)果幾乎沒有差異。以及帶寬的選擇相比較于數(shù)據(jù)的稀少顯得太寬。我們可以定義我們的帶寬函數(shù)以獲得一個更少平滑的結(jié)果。>>>>>>defmy_kde_bandwidth(obj,fac=1./5):..."""WeuseScott'sRule,multipliedbyaconstantfactor."""...returnnp.power(obj.n,-1./(obj.d+4))*fac>>>>>>fig=plt.figure()>>>ax=fig.add_subplot(111)>>>>>>ax.plot(x1,np.zeros(x1.shape),'b+',ms=20)#rugplot>>>kde3=stats.gaussian_kde(x1,bw_method=my_kde_bandwidth)>>>ax.plot(x_eval,kde3(x_eval),'g-',label="WithsmallerBW")>>>>>>plt.show()(Sourcecode)我們看到假如我們設(shè)立帶寬為非常狹窄，則獲得PDF的估計退化為圍繞在數(shù)據(jù)點的簡樸的高斯和。我們現(xiàn)在使用更真實的例子，并且看看在兩種帶寬選擇規(guī)則中的差異。這些規(guī)則被認為在正態(tài)分布上很好用，但即使是偏離正態(tài)的單峰分布上它也工作的很好。作為一個非正態(tài)分布，我們采用5自由度的t分布。importnumpyasnpimportmatplotlib.pyplotaspltfromscipyimportstatsnp.random.seed(12456)x1=np.random.normal(size=200)#randomdata,normaldistributionxs=np.linspace(x1.min()-1,x1.max()+1,200)kde1=stats.gaussian_kde(x1)kde2=stats.gaussian_kde(x1,bw_method='silverman')fig=plt.figure(figsize=(8,6))ax1=fig.add_subplot(211)ax1.plot(x1,np.zeros(x1.shape),'b+',ms=12)#rugplotax1.plot(xs,kde1(xs),'k-',label="Scott'sRule")ax1.plot(xs,kde2(xs),'b-',label="Silverman'sRule")ax1.plot(xs,stats.norm.pdf(xs),'r--',label="TruePDF")ax1.set_xlabel('x')ax1.set_ylabel('Density')ax1.set_title("Normal(top)andStudent'sT$_{df=5}$(bottom)distributions")ax1.legend(loc=1)x2=stats.t.rvs(5,size=200)#randomdata,Tdistributionxs=np.linspace(x2.min()-1,x2.max()+1,200)kde3=stats.gaussian_kde(x2)kde4=stats.gaussian_kde(x2,bw_method='silverman')ax2=fig.add_subplot(212)ax2.plot(x2,np.zeros(x2.shape),'b+',ms=12)#rugplotax2.plot(xs,kde3(xs),'k-',label="Scott'sRule")ax2.plot(xs,kde4(xs),'b-',label="Silverman'sRule")ax2.plot(xs,stats.t.pdf(xs,5),'r--',label="TruePDF")ax2.set_xlabel('x')ax2.set_ylabel('Density')plt.show()(Sourcecode)下面我們看到這個一個寬一個窄的雙峰分布。可以想到結(jié)果將難達成以十分近似，由于每個峰需要不同的帶寬去擬合。>>>>>>fromfunctoolsimportpartial>>>>>>loc1,scal