Statistics(scipy.stats)

Statistics(scipy.stats)................................................................................................................I

介紹..........................................................................2

隨機(jī)變量......................................................................2

獲得幫助.................................................................2

通用方法.................................................................4

位移與縮放...............................................................6

形態(tài)參數(shù).................................................................8

凍結(jié)分布.................................................................9

廣播.....................................................................10

離散分布的特別之處......................................................II

分布擬合.................................................................13

性能問(wèn)題與留意事項(xiàng)......................................................13

遺留問(wèn)題.................................................................14

構(gòu)造詳細(xì)的分布..............................................................14

創(chuàng)建一個(gè)連續(xù)分布，繼承rv_continuous類(lèi)................................14

繼承rv_discrete類(lèi)......................................................16

樣本分析.....................................................................21

描述統(tǒng)計(jì)................................................................21

T檢驗(yàn)和KS檢驗(yàn).........................................................23

分布尾部................................................................25

正態(tài)分布的特別檢驗(yàn)......................................................28

比較兩個(gè)樣本................................................................29

均值.....................................................................30

對(duì)于兩個(gè)不同的樣本進(jìn)行的KS檢驗(yàn)........................................30

核密度估計(jì)..................................................................31

單元估計(jì)................................................................31

多元估計(jì)................................................................40

介紹

在這個(gè)教程我們探討一部分scipy.stats模塊的特性。這里我們的意圖是供應(yīng)應(yīng)運(yùn)用者一個(gè)關(guān)于這個(gè)

包的好用性學(xué)問(wèn)。我們舉薦referencemanual來(lái)介紹更多的細(xì)微環(huán)節(jié)。

留意：這個(gè)文檔還在發(fā)展中<

隨機(jī)變量

有一些通用的概率分布類(lèi)被封裝在continuousrandomvariables以及discreterandomvariables

中。有80多個(gè)連續(xù)性隨機(jī)變量(RVs)以及10余個(gè)離散隨機(jī)變量已經(jīng)用這些類(lèi)建立。同樣，新的程序和

分布可以被用戶(hù)新建(假如你構(gòu)造了一個(gè)，請(qǐng)供應(yīng)它給我們幫助發(fā)展這個(gè)包)。

全部統(tǒng)計(jì)函數(shù)被放在子包scipy.stats中，且有這些函數(shù)的一個(gè)幾乎完整的列表可以運(yùn)用info(stats)

獲得。這個(gè)列表里的隨機(jī)變量也可以從stats子包的docstring中獲得介紹。

在接下來(lái)的探討中，我們著重于連續(xù)性隨機(jī)變量(RVs)。幾乎全部離散變量也符合下面的探討，但是

我們也要指出一些區(qū)分在“離散分布的特別之處”中。

獲得幫助

全部分布可以運(yùn)用help函數(shù)得到說(shuō)明。為獲得這些信息只須要運(yùn)用像這樣的簡(jiǎn)潔調(diào)用：

>>>

>>>fromscipyimportstats

>>>fromscipy.statsimportnorm

>>>printnorm._doc_

作為例子，我們用這種方式找分布的上下界

>>>

>>>print'boundsofdistributionlower:%s,upper:%C%(norm.a.norm.b)

boundsofdistributionlov/er:-infupper:inf

我們可以通過(guò)調(diào)用dir(nonr)來(lái)獲得關(guān)于這個(gè)(正態(tài))分布的全部方法和屬性。應(yīng)當(dāng)看到，一些方法

是私有方法盡管其并沒(méi)有以名稱(chēng)表示出來(lái)(比如它們前面沒(méi)有以卜劃線(xiàn)開(kāi)頭)，比如veccdf就只用于

內(nèi)部計(jì)算(試圖運(yùn)用那些方法將引發(fā)警告，因?yàn)樗鼈兛赡軙?huì)在后續(xù)開(kāi)發(fā)中被移除)

為了獲得真正的主要方法，我們列舉凍結(jié)分布的方法(我們將在下文說(shuō)明何謂“凍結(jié)分布”)

>>>

>>>rv=norm()

>>>dir(rv)#reformatted

['_class_','_delattr_','_dict_'f*_doc_','_getattribute_',

'_hash_','_module_','_new_\'_reduce_'f'_reduce_ex_',

'_repr_','_setattr_','_str_\'_weakref_','args','cdf,'dist',

,entropy','isf,'kwds','moment','pdf','pmf,'ppf','rvs','sf,'stats']

最終，我們能通過(guò)內(nèi)省獲得全部的可用分布的信息。

>>>

>>>norm.cdf(O)

0.5

為了計(jì)算在一個(gè)點(diǎn)上的cdf,我們可以傳遞一個(gè)列表或一個(gè)numpy數(shù)組。

>>>

>>>norm.cdf([-l.,0,1])

array([0.15865525,0.5,0.84134475])

>>>importnumpyasnp

>>>norm.cdf(np.array([-l.,0,1]))

array([0,15865525,0.5,0.841344乃])

相應(yīng)的，像pdf,cdf之類(lèi)的簡(jiǎn)潔方法可以用np.vectorize矢量化.

其他好用的方法可以像這樣運(yùn)用。

>>>

>>>norm.meanO,norm.std(),norm.var()

(0.0,1.0,1.0)

>>>norm.stats(moments="mv")

(array(O.O),array(l.O))

為了找到一個(gè)分布的中心，我們可以運(yùn)用分位數(shù)函數(shù)ppf，它是cdf的逆。

>>>norm.ppf(0.5)

0.0

為了產(chǎn)生一個(gè)隨機(jī)變量集合,

>>>

>>>norm.rvs(size=5)

array([-0.35687759,1.34347647,-0.11710531,-1.00725181,-0.51275702])

不要認(rèn)為norm.rvs(5)產(chǎn)生了五個(gè)變量。

>>>

>>>norm.rvs(5)

欲知其意，請(qǐng)看下一部分的內(nèi)容.

位移與縮放

全部連續(xù)分布可以操縱loc以及scale參數(shù)作為修正location和scale的方式。作為例子，標(biāo)上正態(tài)

分布的location是均值而scale是標(biāo)準(zhǔn)差。

>>>norm.stats(loc=3,scale=4,moments="mv")

(array(3.O),array(16.0))

通常經(jīng)標(biāo)準(zhǔn)化的分布的隨機(jī)變量X可以通過(guò)變換(X-loc)/scale獲得。它們的默認(rèn)值是loc=0以及

scale=l.

聰慧的運(yùn)用loc與scale可以幫助以敏捷的方式調(diào)整標(biāo)準(zhǔn)分布達(dá)到所想要的效果。為了進(jìn)?步說(shuō)明縮

放的效果，下面給出期望為1/人指數(shù)分布的cdf。

F(x)=l-exp(-Xx)

通過(guò)像上面那樣運(yùn)用scale,可以看到得到想要的期望值。

>>>

>>>fromscipy.statsimportexpon

>>>expon.mean(scale=3.)

3.0

勻稱(chēng)分布也是令人感愛(ài)好的：

>>>

>>>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)系起我們?cè)谇懊娑温渲辛粝碌膎orm.rvs(5)的問(wèn)題。事實(shí)上，像這樣調(diào)用?個(gè)分布，其第一

個(gè)參數(shù)，像之前的5,是把loc參數(shù)調(diào)到了5,讓我們看：

>>>

>>>np.mean(norm.rvs(5,size=500))

4.983550784784704

在這里，為說(shuō)明最終一段的輸出：norm.rvs(5)產(chǎn)生了一個(gè)正態(tài)分布變量，其期望，BPloc=5.

我傾向于明確的運(yùn)用loc,scale作為關(guān)鍵字而非像上面那樣利用參數(shù)的依次。因?yàn)檫@看上去有點(diǎn)令人

困惑。我們?cè)谖覀冋f(shuō)明“凍結(jié)RV”的主題之前澄清這一點(diǎn)。

形態(tài)參數(shù)

雖然一個(gè)一般的連續(xù)隨機(jī)變量可以通過(guò)給予loc和scale參數(shù)確定，但一些分布還須要額外的形態(tài)參

數(shù)。作為例子，看到這個(gè)伽馬分布，這是它的密度函數(shù)

y(x,a)=x(Xx)a-ir(a)e-xxz

要求一個(gè)形態(tài)參數(shù)a。留意到X的設(shè)置可以通過(guò)設(shè)置scale關(guān)鍵字為1/人進(jìn)行。

讓我們檢查伽馬分布的形態(tài)參數(shù)的名字的數(shù)量。(我們知道從上面知道其應(yīng)當(dāng)為1)

>>>fromscipy.statsimportgamma

>>>gamma.numargs

1

>>>gamma.shapes

'a'

現(xiàn)在我們?cè)O(shè)置形態(tài)變量的值為1以變成指數(shù)分布。所以我們可以簡(jiǎn)潔的比較是否得到了我們所期望的

結(jié)果。

>>>

>>>gamma(l,scale=2.).stats(moments="mv")

(array(2.0),array(4.0))

留意我們也可以以關(guān)鍵字的方式指定形態(tài)參數(shù)：

>>>

>>>gamma(a=l,scale=2.).stats(moments="mv")

(array(2.0),array(4.0))

凍結(jié)分布

不斷地傳遞loc與scale關(guān)鍵字最終會(huì)讓人厭煩。而凍結(jié)RV的概念被用來(lái)解決這個(gè)問(wèn)題。

>>>rv=gamma(l,scale=2.)

通過(guò)運(yùn)用rv,在任何狀況下我們不再須要包含scale與形態(tài)參數(shù)。明顯，分布可以被多種方式運(yùn)用，

我們可以通過(guò)傳遞全部分布參數(shù)給對(duì)方法的每次調(diào)用(像我們之前做的那樣)或者可以對(duì)一個(gè)分布對(duì)

象凍結(jié)參數(shù)。讓我們看看是懇么回事：

>>>

>>>rv.mean(),rv.std()

(2.0,2.0)

這正是我們應(yīng)當(dāng)?shù)玫降摹?/p>

廣播

像pdf這樣的簡(jiǎn)潔方法滿(mǎn)意numpy的廣播規(guī)則。作為例子，我們可以計(jì)算t分布的右尾分布的臨界

值對(duì)于不同的概率值以及自由度。

>>>

>>>stats.t.isf([0.1,0.05,C.01],[[10],[11]])

array([[1.37218364,1.81246112,2.76376946],

[136343032,1.79588482,2.71807918]])

這里，第一行是以io自由度的臨界值，而其次行是以n為自由度的臨界值。所以，廣播規(guī)則與下面

調(diào)用了兩次isf產(chǎn)生的結(jié)果相同。

>>>stats.t.isf([0.1,0.05,C.01],10)

array([1.37218364,1.81246112,2.76376946])

>>>stats.t.isf([0.1,0.05,C.01],11)

array([1.36343032,1.79588482,2.71807918])

但是假如概率數(shù)組，如[0.1,0.05,0.01]與自由度數(shù)組，如[10,11,12]具有相同的數(shù)組形態(tài)，則進(jìn)行對(duì)應(yīng)

匹配，我們可以分別得到10%,5%,1%尾的臨界值對(duì)于10,11,12的自由度。

>>>

>>>statstisf([O.L0.05,C.01],[10,11,12])

array((1.37218364,1.79538482,2.68099799])

離散分布的特別之處

離散分布的方法的大多數(shù)與連續(xù)分布很類(lèi)似。當(dāng)然像pdf被更換為密度函數(shù)pmf,沒(méi)有估計(jì)方法，

像fit就不能用了。而scale不是一個(gè)合法的關(guān)鍵字參數(shù)。Location參數(shù)，關(guān)鍵字loc則仍舊可以運(yùn)

用用于位移。

cdf的計(jì)算要求一些額外的關(guān)注。在連續(xù)分布的狀況下，累積分布函數(shù)在大多數(shù)標(biāo)準(zhǔn)狀況下是嚴(yán)格遞

增的，所以有唯一的逆。而cdf在離散分布，無(wú)論如何，是階躍函數(shù)，所以cdf的逆，分位點(diǎn)函數(shù)，

要求一個(gè)不同的定義：

ppf(q)=min{x:cdf(x)>=q,xinteger}

為了更多信息口J以看這里.

我們可以看這個(gè)超幾何分布的例子

>>>

>>>fromscipy.statsimporthypergeom

>>>[M,n,N]=[20,7,12]

假如我們?cè)谝恍┱麛?shù)點(diǎn)運(yùn)用cdf,則它們的cdf值再作用ppf會(huì)回到起先的值。

>>>

>>>x=np.arange(4)*2

>>>x

array([0,2,4,6])

>>>prb=hypergeom.cdf(x,M,n,N)

>>>prb

0.9897832817337386])

>>>hypergeom.ppf(prbzM,n,N)

array((0.,2.14.,6.])

假如我們運(yùn)用的值不是cdf的函數(shù)值，則我們得到一個(gè)更高的值。

>>>

>>>hypergeom.ppf(prb-le-8,M,n,N)

array([1.,3,,5.,7.])

>>>hypergeom.ppf(prb-le-8,M,n,N)

array([0.,2.14.,6.])

分布擬合

非凍結(jié)分布的參數(shù)估計(jì)的主要方法：

fit：分布參數(shù)的極大似然估計(jì),包括location與scale

fit.loc.scale:給定形態(tài)參數(shù)確定下的location和scale參數(shù)的估計(jì)

nnlf:負(fù)對(duì)數(shù)似然函數(shù)

expect:計(jì)算函數(shù)pdf或pmf的期望值。

性能問(wèn)題與留意事項(xiàng)

每個(gè)方法的性能與運(yùn)行速度依據(jù)分布的不同表現(xiàn)差異極大。方法的結(jié)果可以由兩種方式獲得，精確計(jì)

和或運(yùn)用獨(dú)立于各詳細(xì)分布的通用算法。

對(duì)于精確計(jì)算，一個(gè)分布能運(yùn)用這種方式的第一種狀況，這個(gè)分布是包中干脆給你的（如標(biāo)準(zhǔn)正態(tài)分

布），其次，給出解析形式，第三通過(guò)scipy.special或numpy.special或numpy.random的rvs特

別函數(shù)給出。一般來(lái)說(shuō)運(yùn)用精確計(jì)算會(huì)比較快。

另一方面，通用方法被用于當(dāng)分布沒(méi)有被指派明確的計(jì)算方法時(shí)運(yùn)用。為了定義一個(gè)分布，只有pdf

或cdf是必需的；通用方法運(yùn)用數(shù)值積分和求根法得到結(jié)果。作為例子，rgh=

stats.gausshyper.rvs（0.5,2,2,2,size=100）以這種方式創(chuàng)建了100個(gè)隨機(jī)變量（抽了100個(gè)值），這

在我的電腦上花了19秒（譯者：我花了3.5秒），對(duì)比取一百萬(wàn)個(gè)標(biāo)準(zhǔn)正態(tài)分布的值只須要1秒。

遺留問(wèn)題

scipy.stats里的分布最近進(jìn)吁了升級(jí)并且被細(xì)致的檢查過(guò)了，不過(guò)仍有一些問(wèn)題存在。

?分布在許多參數(shù)區(qū)間上的值被測(cè)試過(guò)了，然而在一些奇葩的邊界，仍舊可能有錯(cuò)誤的值存在。

?fit的極大似然估計(jì)以默認(rèn)值作為初始參數(shù)將不會(huì)工作的很好，用戶(hù)必需指派合適的初始參

數(shù)。并且，對(duì)于一些分布運(yùn)用極大似然估計(jì)本身就不是一個(gè)好的選擇。

構(gòu)造詳細(xì)的分布

下一個(gè)例子展示了如何建立你自己的分布。更多的例子見(jiàn)分布用法以及統(tǒng)計(jì)檢驗(yàn)

創(chuàng)建一個(gè)連續(xù)分布，繼承rv_continuous類(lèi)

創(chuàng)建連續(xù)分布是特別簡(jiǎn)潔的<

>>>

>>>fromscipyimportstats

>>>classdeterministic_gen(stats.rv_continuous):

...def_cdf(selfx):

...returnnp.where(x<0,0.,1.)

...def_stats(self):

...return0.,0.,0.,0.

>>>

>>>deterministic=dete,ministic_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也能被目動(dòng)計(jì)算出來(lái)：

>>>

>>>deterministic.pdf(np.arange(-3,3,0.5))

array([O.OOOOOOOOe+OO,0.00000000e+00,O.OOOOOOOOe+OO,

O.OOOOOOOOe+00,O.OOOOOOOOe+00,O.OOOOOOOOe+00,

5.83333333e+04,4.16333634e-12,4.16333634e-12,

4.16333634e-12,4.16333634e-12,4.16333634e-12])

留意這種用法的性能問(wèn)題，參見(jiàn)“性能問(wèn)題與留意事項(xiàng)”一節(jié)。這種缺乏信息的通用計(jì)算可能特別慢。

而且再看看?下面這個(gè)精確性的例子：

>>>fromegrateimportquad

>>>quad(deterministic.pdt-le-1,le-1)

但這并不是對(duì)pdf積分的正獺的結(jié)果，事實(shí)上結(jié)果應(yīng)為1.讓我們將積分變得更小一些。

>>>

>>>quad(deterministic.pdf-le-3,le-3)#warningremoved

這樣看上去好多了，然而，問(wèn)題本身來(lái)源于pdf不是來(lái)自包給定的類(lèi)的定義。

繼承rv_discrete類(lèi)

在之后我們運(yùn)用stats.rv_discrete產(chǎn)生一個(gè)離散分布，其有一個(gè)整數(shù)區(qū)間截?cái)喔怕省?/p>

通用信息

通用信息可以從rv_discrete的docstring中得到。

>>>

>>>fromscipy.statsimportrv_discrete

>>>help(rv_discrete)

我們學(xué)到這一點(diǎn)：

“你可以構(gòu)建隨意一個(gè)像P(X=xk)=pk一樣形式的離散rv,通過(guò)傳遞(xk,pk)元組序列給rv_discrete

初始化方法(通過(guò)value=keyword方式)，但其不能有0概率值?！?/p>

接下來(lái)，還有一些進(jìn)一步的要求：

關(guān)鍵字名字必需知道。

Xk點(diǎn)必需是整數(shù)

小數(shù)的有效位數(shù)應(yīng)當(dāng)被給出,

事實(shí)上，假如最終兩個(gè)要求沒(méi)有被滿(mǎn)意，一個(gè)異樣將被拋出或者導(dǎo)致一個(gè)錯(cuò)誤的數(shù)值。

一個(gè)例子

讓我們起先辦，首先

>>>

>>>npoints=20#numberofintegersupportpointsofthedistributionminus1

>>>npointsh=npoints/2

>>>npointsf=float(nponts)

>>>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[:-l]

>>>probs=np.diff(stats.truncnorm.cdf(gridlimitsnorm,-normbound,normbound))

>>>gridint=grid

最終我們可以繼承rv_discrete類(lèi)

>>>

>>>normdiscrete=stats.rv_discrete(vaIues=(gridint,

np.round(probs,decimals=7)),name='normdiscrete')

現(xiàn)在我們已經(jīng)定義了這個(gè)分布，我們可以調(diào)用其全部常規(guī)的離散分布方法。

>>>printmean=%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'))

測(cè)試上面的結(jié)果

讓我們產(chǎn)生一個(gè)隨機(jī)樣本并且比較連續(xù)概率的狀況。

>>>

>>>n_sample=500

>>>np.random.seed(87655678)#fixtheseedforreplicability

>>>rvs=normdiscrete.r^s(size=n_sample)

>>>rvsnd=rvs

>>>fI=np.histogram(r^s,bins=gridlimits)

>>>sfreq=np.vstack([gridint,f,probs*n_sample]).T

>>>printsfreq

[[-1.00000000e+010.0C000000e+002.95019349e-021

[-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+0L]

[-2.00000000e+005.10000000e+015.78004747e+01]

[-1.00000000e+007.10000000e+017.32455414e+01]

[0.00000000e+007.40000000e+017.92618251e+01]

[1.00000000e+008.90000000et017.32455414e+01]

[2.00000000e+005.50000000e+015.78004747e+01]

[3.000000006+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+0D]

[8.00000000e+000.00000000e+005.06497902e-01]

[9.00000000e+000.00000000e+001.32294142e-01]

[1.00000000e+010.00000000e+002.95019349e-02]]

(Sourcecode)

j

u

u

o

n

b

o

上

(Sourcecode)

CumulativeFrequencyandCDFofnormdiscretc

0.6

J

■8

04

0.2

0.0

-10-9-8-7-6-54-3-2-1012345678910

接下來(lái)，我們可以測(cè)試，是否我們的樣本取自于一個(gè)normdiscrete分布。這也是在驗(yàn)證是否陲機(jī)數(shù)

是以F確的方式產(chǎn)生的.

卡方測(cè)試要求至少在每個(gè)子區(qū)間（bin）里具有最小數(shù)目的觀(guān)測(cè)值。我們組合末端子區(qū)間進(jìn)大子區(qū)間

所以它們現(xiàn)在包含了足夠數(shù)量的觀(guān)測(cè)值。

>>>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.4090

P值在這個(gè)狀況下是不顯著地，所以我們可以斷言我們的隨機(jī)樣本的確是由此分布產(chǎn)生的。

樣本分析

首先，我們創(chuàng)建一些隨機(jī)變量。我們?cè)O(shè)置一個(gè)種子所以每次我們都可以得到相同的結(jié)果以便視察。作

為一個(gè)例子，我們從t分布中抽一個(gè)樣本。

>>>

>>>np.random.seed(282629734)

>>>x=stats.t.rvs(10,size=1000)

這里，我們?cè)O(shè)置了t分布的形態(tài)參數(shù)，在這里就是自由度，設(shè)為10.運(yùn)用size=1000表示我們的樣本

由1000個(gè)抽樣是獨(dú)立的(偽。當(dāng)我們不指派loc和scale時(shí)，它們具有默認(rèn)值0和1.

描述統(tǒng)計(jì)

X是一個(gè)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)

如何比較分布本身和它的樣本的指標(biāo)？

>>>

>>>m,v,k=mompnts='mv<;k')

>>>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用的是無(wú)偏的方差估計(jì)量，而np.var卻用的是有偏的估計(jì)量。

T檢驗(yàn)和KS檢驗(yàn)

我們可以運(yùn)用t檢驗(yàn)是否樣本與給定均值(這里是理論均值)存在統(tǒng)計(jì)顯著差異。

>>>

>>>print't-statistic=%6.3fpvalue=%6.4f%stats.ttest_lsamp(x,m)

t-statistic=0.391pvalue=0.6955

P值為0.7,這代表第一類(lèi)錯(cuò)誤的概率，在例子中，為10%。我們不能拒絕“該樣本均值為0”這個(gè)假

設(shè)，0是標(biāo)準(zhǔn)t分布的理論均值。

>>>

>>>tt=(sm-m)/np.sqrt(sv/float(n))#t-statisticformean

>>>pval=stats.t.sf(np.abs(tt),n-l)*2#two-sidedpvafue=Prob(abs(t)>tt)

>>>print't-statistic=%6.3fpvalue=%6.4f%(tt,pvsl)

t-statistic=0391pvalue=0.6955

這里.Kolmogorov-Smirnov檢驗(yàn)(KS檢驗(yàn))被被用來(lái)檢驗(yàn)櫛本是否來(lái)自一個(gè)標(biāo)準(zhǔn)的t分布。

>>>

>>>print'KS-statisticD=%63fpvalue=%6.4f%stats.kstest(x,'t,,(10,))

KS-statisticD=0.016pvalue=0.9606

又?次得到了很高的P值。所以我們不能拒絕樣本是來(lái)自t分布的假設(shè)。在實(shí)際應(yīng)用中，我們不能知

道潛在的分布究竟是什么。假如我們運(yùn)用KS檢驗(yàn)我們的樣本比照正態(tài)分布，那么我們將也不能拒絕

我們的樣本是來(lái)自正態(tài)分布的，在這種狀況下P值為0.4左右,

>>>

>>>print'KS-statisticD=%6.3fpvalue=%6.4f%stats.kstest(x/norm')

KS-statisticD=0.028pvalue=0.3949

無(wú)論如何，標(biāo)準(zhǔn)正態(tài)分布行1的方差，當(dāng)我們的樣本有1.29時(shí)，假如我們標(biāo)準(zhǔn)化我們的樣本并且測(cè)試

它比照正態(tài)分布，那么P值將又一次很高我們將還是不能拒絕假設(shè)是來(lái)自正態(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

注釋?zhuān)篕S檢驗(yàn)假設(shè)我們比照的分布就是以給定的參數(shù)確定的，但我們?cè)谧罱K估計(jì)了均值和方差，這

個(gè)假設(shè)就被違反了，故而這個(gè)測(cè)試統(tǒng)計(jì)量的P值是含偏的，這個(gè)用法是錯(cuò)誤的。

分布尾部

最終，我們可以檢查分布的石尾，我們可以運(yùn)用分位點(diǎn)函數(shù)ppf，其為cdf函數(shù)的逆，來(lái)獲得臨界值，

或者更干脆的，我們可以運(yùn)用殘存函數(shù)的逆來(lái)辦。

>>>

>>>critOl,critOS,critlO=stats.t.ppf([l-O.Ol,1-0.05,1-0.10],10)

>>>print'criticalvaluesfromppfat1%%,5%%and10%%%8.4f%8.4f%8.4f%(critOl,c

ritO5,critlO)

criticalvaluesfromppfat1%,5%and10%2.76381.81251.3722

>>>print'criticalvaluesfromisfat1%%,5%%and10%%%8.4f%8.4f%8.4f%tuple(stat

s.t.isf([0.01,0.05,0.10],10))

criticalvaluesfromisfat1%,5%and10%2.76381.81251.3722

>>>

>>>freqOl=np.sum(x>crit01)/float(n)*100

>>>freqOS=np.sum(x>crit05)/float(n)*100

>>>freq10=np.sum(x>critl0)/float(n)*100

>>>print'sample%^-frequencyat1%%,5%%and10%%tail%8.4f%8.4f%8.4f%(freqO

1,freqOS,freq10)

sample%-frequencyat1%,5%and10%tail1.40005.800010.5000

在這三種狀況中，我們的樣本有有一個(gè)更重的尾部，即實(shí)際在理論分界值右邊的概率要高于理論值。

我們可以通過(guò)運(yùn)用更大的樣原來(lái)獲得更好的擬合。在以下?tīng)顩r閱歷頻率已經(jīng)很接近理論概率了，但即

使我們重復(fù)這個(gè)過(guò)程若干次，波動(dòng)照舊會(huì)保持在這個(gè)程度。

>>>freqO5l=np.sum(stats.t.rvs(10,size=10000)>critOS)/10000.0*100

>>>print'largersample^%-frequencyat5%%tail%8.4f%freq051

largersample%-frequencyat5%tail4.8000

我們也可以比較它與正態(tài)分布的尾部，其有一個(gè)輕的多的尾部：

>>>

>>>print'tailprob,ofnormalat1%%,5%%and10%%%8.4f%8.4f%8.4f%\

...norm?;f([rritO1,rritOS,rrit1O])*1OO)

tailprob,ofnormalat1%,5%and10%0.28573.49578.5003

卡方檢驗(yàn)可以被用來(lái)測(cè)試，是否?個(gè)有限的分類(lèi)觀(guān)測(cè)值頻率與假定的理論概率分布具有顯著差異。

>>>

>>>quantiles=[0.0,0.01,0.05,0.1,1-0.10,1-0.05,1-0.01,1.0]

>>>crit=stats.t.ppf(quaitiles,10)

>>>printcrit

[-Inf-2.76376946-1.81246112-1.372183641.372183641.81246112

2.76376946Inf]

>>>n_sample=x.size

>>>freqcount=np.histogram(x,bins=crit)[O]

>>>tprob=np.diff(quantiles)

>>>nprob=np.diff(stats.norm.cdf(crit))

>>>tch,tpval=stats.chisquare(freqcounttprob*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分相檢驗(yàn)沒(méi)被拒絕時(shí)標(biāo)準(zhǔn)正態(tài)分布卻被完全拒絕。在我們的樣本區(qū)分出這兩個(gè)分布后，

我們可以先進(jìn)行擬合確定scale與location再檢杳擬合后的分布的差異性。

我們可以先進(jìn)行擬合，再用擬合分布而不是默認(rèn)(至少location和scale是默認(rèn)的)分布去進(jìn)行檢驗(yàn)。

>>>

>>>tdoftloc,tscale=stats.t.fit(x)

>>>nloc,nscale=stats.norm.fit(x)

>>>tprob=np.diff(stats.t.cdf(crit,tdofloc=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

在經(jīng)過(guò)參數(shù)調(diào)整之后，我們?nèi)耘f可以以5%水平拒絕正態(tài)分布假設(shè)。然而卻以95%的p值明顯的不能

拒絕t分布。

正態(tài)分布的特別檢驗(yàn)

自從正態(tài)分布變?yōu)榻y(tǒng)計(jì)學(xué)中最常見(jiàn)的分布，就出現(xiàn)了大量的方法用來(lái)檢驗(yàn)一個(gè)樣本是否可以被看成是

來(lái)自正態(tài)分布的。

苜先我們檢驗(yàn)分布的峰度和倍度是否顯著地與正態(tài)分布的對(duì)應(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

將這兩個(gè)檢驗(yàn)組合起來(lái)的正態(tài)性檢驗(yàn)

>>>

>>>print'normaltestteststat=%6.3fpvalue=%6.4f%stats.normaltest(x)

normaltestteststat=30.379pvalue=0.0000

在全部三個(gè)測(cè)試中，P值是特別低的，所以我們可以拒絕我們的樣本的峰度與偏度與正態(tài)分布相同的

假設(shè)。

當(dāng)我們的樣本標(biāo)準(zhǔn)化之后，我們照舊得到相同的結(jié)果。

>>>

>>>print'normaltestteststat=%6.3fpvalue=%6.4f%\

...stats.normaltest((x-x.mean())/x.stdi))

normaltestteststat=30.379pvalue=0.0000

因?yàn)镕態(tài)件被很強(qiáng)的拒絕了,所以我們可以檢查這種檢驗(yàn)方式是否可以有效地作用到其他狀況中°

>>>

>>>print'normaltestteststat=%6.3fpvalue=%6.4f%stats.normaltest(stats.t.rvs(13,siz

e=100))

normaltestteststat=4.698pvalue=0.0955

>>>print'normaltestteststat=%6.3fpvalue=%6.4f%stats.normaltest(stats.norm.r^s(siz

e=1000))

normaltestteststat=0.613pvalue=0.7361

我們檢驗(yàn)了小樣本的t分布樣本的觀(guān)測(cè)值以及一個(gè)大樣本的正態(tài)分布觀(guān)測(cè)值，在這兩種狀況中我們都

不能拒絕其來(lái)自正態(tài)分布的空假設(shè)。得到這樣的結(jié)果是因?yàn)榍罢呤且驗(yàn)闊o(wú)法區(qū)分小樣本下的t分布，

后者是因?yàn)樗瓉?lái)就來(lái)自正態(tài)分布。

比較兩個(gè)樣本

接下來(lái)，我們有兩個(gè)分布，其可以判定為相同或者來(lái)自不同的分布，以及我們希望測(cè)試是否這些樣本

有相同的統(tǒng)計(jì)特征。

均值

以相同的均值產(chǎn)生的樣本進(jìn)行檢驗(yàn)：

>>>

>>>rvsl=stats.norm.rvs(loc=5,scale=10,size=500)

>>>rvs2=stats.norm.rvs(loc=5,scale=10,size=500)

>>>stats.ttest_ind(rvsl,rvs2)

以不同的均值產(chǎn)生的樣本進(jìn)行檢驗(yàn)：

>>>

>>>rvs3=stats.norm.rvs(loc=8,scale=10,size=500)

>>>stats.ttest_ind(rvsl,rvs3)

對(duì)于兩個(gè)不同的樣本進(jìn)行的KS檢驗(yàn)

在這個(gè)例子中我們運(yùn)用兩個(gè)同分布的樣本進(jìn)行檢驗(yàn).設(shè)因?yàn)镻值很高，亳不驚奇我們不能拒絕原假。

>>>stats.ks_2samp(rvsl,rvs2)

在其次個(gè)例子中，由于均值不同，所以我們可以拒絕空假設(shè)，由P值小于1%。

>>>

>>>stats.ks_2samp(rvsl,rvs3)

核密度估計(jì)

一個(gè)常見(jiàn)的統(tǒng)計(jì)學(xué)問(wèn)題是從一個(gè)樣本中估計(jì)隨機(jī)變量的概率密度分布函數(shù)（PDF）這個(gè)問(wèn)題被稱(chēng)為密

度估計(jì)，對(duì)此最聞名的工具是直方圖。直方圖是?個(gè)很好的可視化工具（主要是因?yàn)槊總€(gè)人都理解它）。

但是對(duì)于對(duì)于數(shù)據(jù)特征的利用卻并不是特別有效率。核密度估計(jì)（KDE對(duì)于這個(gè)問(wèn)題）是一個(gè)更有效的

工具。這個(gè)gaussian_kde缶計(jì)方法可以被用來(lái)估計(jì)單元或多元數(shù)據(jù)的PDF。它在數(shù)據(jù)呈單峰的時(shí)候

工作的最好，但也可以在多峰狀況下工作。

單元估計(jì)

我們以一個(gè)最小數(shù)據(jù)集來(lái)視察gaussian_kde是如何工作的，以及帶寬（bandwidth）的不同選擇方

式。PDF對(duì)應(yīng)的數(shù)據(jù)被以藍(lán)線(xiàn)的形式顯示在圖像的底端（被稱(chēng)為毯圖（rugplot））

>>>

>>>fromscipyimportstats

>>>importmatplotlib.pyplotaspit

>>>xl=np.array([-7,-5,1,4,5],dtype=np.float)

>>>kdel=stats.gaussian_kde(xl)

>>>kde2=stats.gaussian_kde(xl,bw_method='silverman')

>>>

>>>fig=plt.figureO

>>>ax=fig.add_subplot(lll)

>>>

>>>ax.plot(xl,np.zeros(xl.shape),'b+',ms=20)#rugplot

>>>x_eval=np.linspace(-10,10,num=200)

>>>ax.plot(x_eval,kdel(x_eval),'k-',label="Scott'sRule")

>>>ax.plot(x_eval,kdel(x_eval),'r-\label="Silverman'sRule")

>>>

>>>plt.showQ

(Sourcecode)

我們看到在Scott規(guī)則以及Silverman規(guī)則下的結(jié)果幾乎沒(méi)有差異。以及帶寬的選擇相比較于數(shù)據(jù)的

稀有顯得太寬。我們可以定義我們的帶寬函數(shù)以獲得一個(gè)更少平滑的結(jié)果。

>>>

>>>defmy_kde_bandwidth(obj,fac=l./5):

...'**144?useScott'sRule,multipliedbyaconstantfactor：""

...returnnp.power(obj.n,-l./(obj.d+4))*fac

>>>

>>>fig=plt.figureO

>>>ax=fig.add_subplot(lll)

>>>

>>>ax.plot(xl,np.zeros(xl.shape),'b+',ms=20)#rugplot

>>>kde3=stats.gaussian_kde(xl,bw_method=my_kde_bandwidth)

>>>ax.plot(x_eval,kde3(x_eval),'g-',label："WithsmallerBW")

>>>

>>>plt.showQ

(Sourcecode)

我們看到假如我們?cè)O(shè)置帶寬為特別狹窄，則獲得PDF的估計(jì)退化為圍繞在數(shù)據(jù)點(diǎn)的簡(jiǎn)潔的高斯和。

我們現(xiàn)在運(yùn)用更真實(shí)的例子，并且看看在兩種帶寬選擇規(guī)則中的差異。這些規(guī)則被認(rèn)為在正態(tài)分布上

很好用，但即使是偏離正態(tài)的單峰分布上它也工作的很好。作為一個(gè)非正態(tài)分布，我們采納5自由度

的t分布。

importnumpyasnp

importmatplotlib.pyplotaspit

fromscipyimportstats

np.random.seed(12456)

xl=np.random.normal(size=200)#randomdata,normaldistribution

xs=np.linspace(xl.min()-l,xl.max()+l,200)

kdel=stats.gaussian_kde[xl)

kde2=stats.gaussian_kde[xl,bw_method='silverman')

fig=plt.figure(figsize=(8,6))

axl=fig.add_subplot(211)

axl.plot(xl,np.zeros(xl.shape),'b+',ms=12)#rugplot

axl.plot(xs,kdel(xs),'k-',label="Scott'sRule")

axl.plot(xs,kde2(xs),'b-',abel="Silverman'sRule')

axl.plot(xs,stats.norm.pdf(xs),'r—label="TruePDF")

axl.set_xlabel('x')

axl.set_ylabel('Density')

axl.set_title("Normal(top)andStudent'sT$_{df=5}$(bottom)distributions")

axl.legend(loc=l)

x2=stats.t.rvs(5,size=200)#randomdata,Tdistribution

xs=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)#rugplot

ax2.plot(xs,kde3(xs),'k-'flabel="Scott'sRule")

ax2.plot(xs,kde4(xs),'b-\abel="Silverman'sRule')

ax2.plot(xs,stats.t.pdf(xs,5),'r—label="TruePDF")

ax2.set_xlabel('x')

ax2.set_ylabel('Density')

plt.showQ

(Sourcecode)

Normal(top)andStudcnl'sT.,>s(bottom)distributions

0.40=

035

030

>

20.25

30.20

0.15

0.10

0.05

0.00

040

035

0.30

0.25

0.20

0.15

0.10

0.05

aoo

下面我們看到這個(gè)一個(gè)寬一個(gè)窄的雙峰分布。可以想到結(jié)果將難達(dá)到以特別近似，因?yàn)槊總€(gè)峰須要不

同的帶寬去擬合。

>>>

>>>fromfunctoolsimportpartial

>>>

>>>loci,