Python統(tǒng)計(jì)學(xué)包scipystats手冊

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

2、提供給使用者一個(gè)關(guān)于這個(gè)包的實(shí)用性知識。我們推薦reference manual來介紹更多的細(xì)節(jié)。注意：這個(gè)文檔還在發(fā)展中。隨機(jī)變量有一些通用的概率分布類被封裝在continuous random variables以及discrete random variables中。有80多個(gè)連續(xù)性隨機(jī)變量(RVs)以及10余個(gè)離散隨機(jī)變量已經(jīng)用這些類建立。同樣，新的程序和分布可以被用戶新建（如果你構(gòu)造了一個(gè)，請?zhí)峁┧o我們幫助發(fā)展這個(gè)包）。所有統(tǒng)計(jì)函數(shù)被放在子包scipy.stats中，且有這些函數(shù)的一個(gè)幾乎完整的列表可以使用info(stats)獲得。這個(gè)列表里的隨機(jī)變量也可以從stats子包的do

3、cstring中獲得介紹。在接下來的討論中，我們著重于連續(xù)性隨機(jī)變量(RVs)。幾乎所有離散變量也符合下面的討論，但是我們也要指出一些區(qū)別在“離散分布的特殊之處”中。獲得幫助所有分布可以使用help函數(shù)得到解釋。為獲得這些信息只需要使用像這樣的簡單調(diào)用： from scipy import stats from scipy.stats import norm print norm._doc_作為例子，我們用這種方式找分布的上下界 print bounds of distribution lower: %s, upper: %s % (norm.a,norm.b)bounds of distri

4、bution lower: -inf, upper: inf我們可以通過調(diào)用dir(norm)來獲得關(guān)于這個(gè)（正態(tài)）分布的所有方法和屬性。應(yīng)該看到，一些方法是私有方法盡管其并沒有以名稱表示出來（比如它們前面沒有以下劃線開頭），比如veccdf就只用于內(nèi)部計(jì)算（試圖使用那些方法將引發(fā)警告，因?yàn)樗鼈兛赡軙诤罄m(xù)開發(fā)中被移除）為了獲得真正的主要方法，我們列舉凍結(jié)分布的方法（我們將在下文解釋何謂“凍結(jié)分布”） rv = norm() dir(rv) # reformatted _class_, _delattr_, _dict_, _doc_, _getattribute_, _hash_, _ini

5、t_, _module_, _new_, _reduce_, _reduce_ex_, _repr_, _setattr_, _str_, _weakref_, args, cdf, dist, entropy, isf, kwds, moment, pdf, pmf, ppf, rvs, sf, stats最后，我們能通過內(nèi)省獲得所有的可用分布的信息。 import warnings warnings.simplefilter(ignore, DeprecationWarning) dist_continu = d for d in dir(stats) if. isinstance(get

6、attr(stats,d), stats.rv_continuous) dist_discrete = d for d in dir(stats) if. isinstance(getattr(stats,d), stats.rv_discrete) print number of continuous distributions:, len(dist_continu)number of continuous distributions: 84 print number of discrete distributions: , len(dist_discrete)number of discr

7、ete distributions: 12通用方法連續(xù)隨機(jī)變量的主要公共方法如下： rvs:隨機(jī)變量（就是從這個(gè)分布中抽一些樣本） pdf：概率密度函數(shù)。 cdf：累計(jì)分布函數(shù) sf：殘存函數(shù)（1-CDF） ppf：分位點(diǎn)函數(shù)（CDF的逆） isf：逆殘存函數(shù)（sf的逆） stats:返回均值，方差，（費(fèi)舍爾）偏態(tài)，（費(fèi)舍爾）峰度。 moment:分布的非中心矩。讓我們使用一個(gè)標(biāo)準(zhǔn)的RV作為例子。 norm.cdf(0)0.5為了計(jì)算在一個(gè)點(diǎn)上的cdf，我們可以傳遞一個(gè)列表或一個(gè)numpy數(shù)組。 norm.cdf(-1., 0, 1)array( 0.15865525, 0.5 , 0.841

8、34475) import numpy as np norm.cdf(np.array(-1., 0, 1)array( 0.15865525, 0.5 , 0.84134475)相應(yīng)的，像pdf,cdf之類的簡單方法可以用np.vectorize矢量化.其他實(shí)用的方法可以像這樣使用。 norm.mean(), norm.std(), norm.var()(0.0, 1.0, 1.0) norm.stats(moments = mv)(array(0.0), array(1.0)為了找到一個(gè)分布的中心，我們可以使用分位數(shù)函數(shù)ppf，它是cdf的逆。 norm.ppf(0.5)0.0為了產(chǎn)生一個(gè)

9、隨機(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)欲知其意，請看下一部分的內(nèi)容。位移與縮放位移與所有連續(xù)分布可以操縱loc以及scale參數(shù)作為修正location和scale的方式。作為例子，標(biāo)準(zhǔn)正態(tài)分布的location是均值而scale是標(biāo)準(zhǔn)差。 norm.stats(loc = 3, scale = 4, moments = mv)(array(3.0), array(16.0)通常

10、經(jīng)標(biāo)準(zhǔn)化的分布的隨機(jī)變量X可以通過變換(X-loc)/scale獲得。它們的默認(rèn)值是loc=0以及scale=1.聰明的使用loc與scale可以幫助以靈活的方式調(diào)整標(biāo)準(zhǔn)分布達(dá)到所想要的效果。為了進(jìn)一步說明縮放的效果，下面給出期望為1/指數(shù)分布的cdf。F(x)=1exp(x)通過像上面那樣使用scale，可以看到得到想要的期望值。 from scipy.stats import expon expon.mean(scale=3.)3.0均勻分布也是令人感興趣的： from scipy.stats import uniform uniform.cdf(0, 1, 2, 3, 4, 5, loc

11、 = 1, scale = 4)array( 0. , 0. , 0.25, 0.5 , 0.75, 1. )最后，聯(lián)系起我們在前面段落中留下的norm.rvs(5)的問題。事實(shí)上，像這樣調(diào)用一個(gè)分布，其第一個(gè)參數(shù)，像之前的5，是把loc參數(shù)調(diào)到了5，讓我們看： np.mean(norm.rvs(5, size=500)4.983550784784704在這里，為解釋最后一段的輸出：norm.rvs(5)產(chǎn)生了一個(gè)正態(tài)分布變量，其期望，即loc=5.我傾向于明確的使用loc,scale作為關(guān)鍵字而非像上面那樣利用參數(shù)的順序。因?yàn)檫@看上去有點(diǎn)令人困惑。我們在我們解釋“凍結(jié)RV”的主題之前澄清這一

12、點(diǎn)。形態(tài)參數(shù)雖然一個(gè)一般的連續(xù)隨機(jī)變量可以通過賦予loc和scale參數(shù)確定，但一些分布還需要額外的形態(tài)參數(shù)。作為例子，看到這個(gè)伽馬分布，這是它的密度函數(shù)(x,a)=(x)a1(a)ex,要求一個(gè)形態(tài)參數(shù)a。注意到的設(shè)置可以通過設(shè)置scale關(guān)鍵字為1/進(jìn)行。讓我們檢查伽馬分布的形態(tài)參數(shù)的名字的數(shù)量。（我們知道從上面知道其應(yīng)該為1） from scipy.stats import gamma gamma.numargs1 gamma.shapesa現(xiàn)在我們設(shè)置形態(tài)變量的值為1以變成指數(shù)分布。所以我們可以容易的比較是否得到了我們所期望的結(jié)果。 gamma(1, scale=2.).stats(m

13、oments=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的概念被用來解決這個(gè)問題。 rv = gamma(1, scale=2.)通過使用rv，在任何情況下我們不再需要包含scale與形態(tài)參數(shù)。顯然，分布可以被多種方式使用，我們可以通過傳遞所有分布參數(shù)給對方法的每次調(diào)用（像我們之前做的那樣）或者可以對一個(gè)分布對象凍結(jié)參數(shù)。讓我們看看是怎么回事： rv.m

14、ean(), rv.std()(2.0, 2.0)這正是我們應(yīng)該得到的。廣播像pdf這樣的簡單方法滿足numpy的廣播規(guī)則。作為例子，我們可以計(jì)算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

15、.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)，則進(jìn)行對應(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.6809979

16、9)離散分布的特殊之處離散分布的方法的大多數(shù)與連續(xù)分布很類似。當(dāng)然像pdf被更換為密度函數(shù)pmf，沒有估計(jì)方法，像fit就不能用了。而scale不是一個(gè)合法的關(guān)鍵字參數(shù)。Location參數(shù)，關(guān)鍵字loc則仍然可以使用用于位移。cdf的計(jì)算要求一些額外的關(guān)注。在連續(xù)分布的情況下，累積分布函數(shù)在大多數(shù)標(biāo)準(zhǔn)情況下是嚴(yán)格遞增的，所以有唯一的逆。而cdf在離散分布，無論如何，是階躍函數(shù)，所以cdf的逆，分位點(diǎn)函數(shù)，要求一個(gè)不同的定義：ppf(q) = minx : cdf(x) = q, x integer為了更多信息可以看這里。我們可以看這個(gè)超幾何分布的例子 from scipy.stats imp

17、ort hypergeom M, n, N = 20, 7, 12如果我們在一些整數(shù)點(diǎn)使用cdf，則它們的cdf值再作用ppf會回到開始的值。 x = np.arange(4)*2 xarray(0, 2, 4, 6) prb = hypergeom.cdf(x, M, n, N) prb 0.9897832817337386) hypergeom.ppf(prb, M, n, N)array( 0., 2., 4., 6.)如果我們使用的值不是cdf的函數(shù)值，則我們得到一個(gè)更高的值。 hypergeom.ppf(prb + 1e-8, M, n, N)array( 1., 3., 5., 7

18、.) hypergeom.ppf(prb - 1e-8, M, n, N)array( 0., 2., 4., 6.)分布擬合非凍結(jié)分布的參數(shù)估計(jì)的主要方法：fit：分布參數(shù)的極大似然估計(jì)，包括location與scalefit_loc_scale: 給定形態(tài)參數(shù)確定下的location和scale參數(shù)的估計(jì)nnlf:負(fù)對數(shù)似然函數(shù)expect:計(jì)算函數(shù)pdf或pmf的期望值。性能問題與注意事項(xiàng)每個(gè)方法的性能與運(yùn)行速度根據(jù)分布的不同表現(xiàn)差異極大。方法的結(jié)果可以由兩種方式獲得，精確計(jì)算或使用獨(dú)立于各具體分布的通用算法。對于精確計(jì)算，一個(gè)分布能使用這種方式的第一種情況，這個(gè)分布是包中直接給你的（如

19、標(biāo)準(zhǔn)正態(tài)分布），第二，給出解析形式，第三通過scipy.special或numpy.special或numpy.random的rvs特殊函數(shù)給出。一般來說使用精確計(jì)算會比較快。另一方面，通用方法被用于當(dāng)分布沒有被指派明確的計(jì)算方法時(shí)使用。為了定義一個(gè)分布，只有pdf或cdf是必須的；通用方法使用數(shù)值積分和求根法得到結(jié)果。作為例子，rgh = stats.gausshyper.rvs(0.5, 2, 2, 2, size=100)以這種方式創(chuàng)建了100個(gè)隨機(jī)變量（抽了100個(gè)值），這在我的電腦上花了19秒（譯者：我花了3.5秒），對比取一百萬個(gè)標(biāo)準(zhǔn)正態(tài)分布的值只需要1秒。遺留問題scipy.st

20、ats里的分布最近進(jìn)行了升級并且被仔細(xì)的檢查過了，不過仍有一些問題存在。 分布在很多參數(shù)區(qū)間上的值被測試過了，然而在一些奇葩的邊界，仍然可能有錯(cuò)誤的值存在。 fit的極大似然估計(jì)以默認(rèn)值作為初始參數(shù)將不會工作的很好，用戶必須指派合適的初始參數(shù)。并且，對于一些分布使用極大似然估計(jì)本身就不是一個(gè)好的選擇。構(gòu)造具體的分布下一個(gè)例子展示了如何建立你自己的分布。更多的例子見分布用法以及統(tǒng)計(jì)檢驗(yàn)創(chuàng)建一個(gè)連續(xù)分布，繼承rv_continuous類創(chuàng)建連續(xù)分布是非常簡單的。 from scipy import stats class deterministic_gen(stats.rv_continuous)

21、:. def _cdf(self, x):. return np.where(x 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也能被自動(dòng)計(jì)算出來： deterministic.pdf(np.arange(-3, 3, 0.5)array( 0.00000000e+00, 0.00000000e+00, 0.00000000e+0

22、0, 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)注意這種用法的性能問題，參見“性能問題與注意事項(xiàng)”一節(jié)。這種缺乏信息的通用計(jì)算可能非常慢。而且再看看下面這個(gè)準(zhǔn)確性的例子： from egrate import quad quad(deterministic.pdf, -1e-1, 1e-1)但這并不是對pdf積分的正確的結(jié)果，實(shí)際上結(jié)果

23、應(yīng)為1.讓我們將積分變得更小一些。 quad(deterministic.pdf, -1e-3, 1e-3) # warning removed(1.000這樣看上去好多了，然而，問題本身來源于pdf不是來自包給定的類的定義。繼承rv_discrete類在之后我們使用stats.rv_discrete產(chǎn)生一個(gè)離散分布，其有一個(gè)整數(shù)區(qū)間截?cái)喔怕?。通用信息通用信息可以?rv_discrete的 docstring中得到。 from scipy.stats import rv_discrete help(rv_discrete)我們學(xué)到這一點(diǎn)：“你可以構(gòu)建任意一個(gè)像P(X=xk)=pk一樣形式的離

24、散rv，通過傳遞(xk,pk)元組序列給rv_discrete初始化方法(通過value=keyword方式)，但其不能有0概率值?！苯酉聛?，還有一些進(jìn)一步的要求：關(guān)鍵字名字必須知道。Xk點(diǎn)必須是整數(shù)小數(shù)的有效位數(shù)應(yīng)當(dāng)被給出。事實(shí)上，如果最后兩個(gè)要求沒有被滿足，一個(gè)異常將被拋出或者導(dǎo)致一個(gè)錯(cuò)誤的數(shù)值。一個(gè)例子讓我們開始辦，首先 npoints = 20 # number of integer support points of the distribution minus 1 npointsh = npoints / 2 npointsf = float(npoints) nbound = 4

25、 # bounds for the truncated normal normbound = (1+1/npointsf) * nbound # actual bounds of truncated normal grid = np.arange(-npointsh, npointsh+2, 1) # integer grid gridlimitsnorm = (grid-0.5) / npointsh * nbound # bin limits for the truncnorm gridlimits = grid - 0.5 # used later in the analysis gri

26、d = 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)定義了這個(gè)分布，我們可以調(diào)用其所有常規(guī)的離散分布方法。 print mean = %6.4f, variance = %6.4f, skew = %

27、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)生一個(gè)隨機(jī)樣本并且比較連續(xù)概率的情況。 n_sample = 500 np.random.seed(87655678) # fix the seed for replicability rvs = normdiscrete.rv

28、s(size=n_sample) rvsnd = rvs f, l = np.histogram(rvs, bins=gridlimits) sfreq = np.vstack(gridint, f, probs*n_sample).T print sfreq -1.00000000e+01 0.00000000e+00 2.95019349e-02 -9.00000000e+00 0.00000000e+00 1.32294142e-01 -8.00000000e+00 0.00000000e+00 5.06497902e-01 -7.00000000e+00 2.00000000e+00

29、1.65568919e+00 -6.00000000e+00 1.00000000e+00 4.62125309e+00 -5.00000000e+00 9.00000000e+00 1.10137298e+01 -4.00000000e+00 2.60000000e+01 2.24137683e+01 -3.00000000e+00 3.70000000e+01 3.89503370e+01 -2.00000000e+00 5.10000000e+01 5.78004747e+01 -1.00000000e+00 7.10000000e+01 7.32455414e+01 0.0000000

30、0e+00 7.40000000e+01 7.92618251e+01 1.00000000e+00 8.90000000e+01 7.32455414e+01 2.00000000e+00 5.50000000e+01 5.78004747e+01 3.00000000e+00 5.00000000e+01 3.89503370e+01 4.00000000e+00 1.70000000e+01 2.24137683e+01 5.00000000e+00 1.10000000e+01 1.10137298e+01 6.00000000e+00 4.00000000e+00 4.6212530

31、9e+00 7.00000000e+00 3.00000000e+00 1.65568919e+00 8.00000000e+00 0.00000000e+00 5.06497902e-01 9.00000000e+00 0.00000000e+00 1.32294142e-01 1.00000000e+01 0.00000000e+00 2.95019349e-02(Source code)(Source code)接下來，我們可以測試，是否我們的樣本取自于一個(gè)normdiscrete分布。這也是在驗(yàn)證是否隨機(jī)數(shù)是以正確的方式產(chǎn)生的?？ǚ綔y試要求起碼在每個(gè)子區(qū)間（bin）里具有最小數(shù)目的觀測

32、值。我們組合末端子區(qū)間進(jìn)大子區(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 chisquare for normdiscrete: chi2 = %6.3f pvalue = %6.4f % (ch2, pval)chisquare for normdiscrete: chi2 = 12.466 pva

33、lue = 0.4090P值在這個(gè)情況下是不顯著地，所以我們可以斷言我們的隨機(jī)樣本的確是由此分布產(chǎn)生的。樣本分析樣本分析yangya首先，我們創(chuàng)建一些隨機(jī)變量。我們設(shè)置一個(gè)種子所以每次我們都可以得到相同的結(jié)果以便觀察。作為一個(gè)例子，我們從t分布中抽一個(gè)樣本。 np.random.seed(282629734) x = stats.t.rvs(10, size=1000)這里，我們設(shè)置了t分布的形態(tài)參數(shù)，在這里就是自由度，設(shè)為10.使用size=1000表示我們的樣本由1000個(gè)抽樣是獨(dú)立的（偽）。當(dāng)我們不指派loc和scale時(shí)，它們具有默認(rèn)值0和1.描述統(tǒng)計(jì)X是一個(gè)numpy數(shù)組。我們可以直

34、接調(diào)用它所有的方法。 print x.max(), x.min() # equivalent to np.max(x), np.min(x)5.26327732981 -3.78975572422 print x.mean(), x.var() # equivalent to np.mean(x), np.var(x)如何比較分布本身和它的樣本的指標(biāo)？ m, v, s, k = stats.t.stats(10, moments=mvsk) n, (smin, smax), sm, sv, ss, sk = stats.describe(x) print distribution:,distr

35、ibution: sstr = mean = %6.4f, variance = %6.4f, skew = %6.4f, kurtosis = %6.4f print sstr %(m, v, s ,k)mean = 0.0000, variance = 1.2500, skew = 0.0000, kurtosis = 1.0000 print sample: ,sample: print sstr %(sm, sv, ss, sk)mean = 0.0141, variance = 1.2903, skew = 0.2165, kurtosis = 1.0556注意：stats.desc

36、ribe用的是無偏的方差估計(jì)量，而np.var卻用的是有偏的估計(jì)量。T檢驗(yàn)和KS檢驗(yàn)我們可以使用t檢驗(yàn)是否樣本與給定均值（這里是理論均值）存在統(tǒng)計(jì)顯著差異。 print t-statistic = %6.3f pvalue = %6.4f % stats.ttest_1samp(x, m)t-statistic = 0.391 pvalue = 0.6955P值為0.7，這代表第一類錯(cuò)誤的概率，在例子中，為10%。我們不能拒絕“該樣本均值為0”這個(gè)假設(shè)，0是標(biāo)準(zhǔn)t分布的理論均值。 tt = (sm-m)/np.sqrt(sv/float(n) # t-statistic for mean pv

37、al = stats.t.sf(np.abs(tt), n-1)*2 # two-sided pvalue = Prob(abs(t)tt) print t-statistic = %6.3f pvalue = %6.4f % (tt, pval)t-statistic = 0.391 pvalue = 0.6955這里Kolmogorov-Smirnov檢驗(yàn)（KS檢驗(yàn)）被被用來檢驗(yàn)樣本是否來自一個(gè)標(biāo)準(zhǔn)的t分布。 print KS-statistic D = %6.3f pvalue = %6.4f % stats.kstest(x, t, (10,)KS-statistic D = 0.01

38、6 pvalue = 0.9606又一次得到了很高的P值。所以我們不能拒絕樣本是來自t分布的假設(shè)。在實(shí)際應(yīng)用中，我們不能知道潛在的分布到底是什么。如果我們使用KS檢驗(yàn)我們的樣本對照正態(tài)分布，那么我們將也不能拒絕我們的樣本是來自正態(tài)分布的，在這種情況下P值為0.4左右。 print KS-statistic D = %6.3f pvalue = %6.4f % stats.kstest(x,norm)KS-statistic D = 0.028 pvalue = 0.3949無論如何，標(biāo)準(zhǔn)正態(tài)分布有1的方差，當(dāng)我們的樣本有1.29時(shí)。如果我們標(biāo)準(zhǔn)化我們的樣本并且測試它比照正態(tài)分布，那么P值將又一

39、次很高我們將還是不能拒絕假設(shè)是來自正態(tài)分布的。 d, pval = stats.kstest(x-x.mean()/x.std(), norm) print KS-statistic D = %6.3f pvalue = %6.4f % (d, pval)KS-statistic D = 0.032 pvalue = 0.2402注釋：KS檢驗(yàn)假設(shè)我們比照的分布就是以給定的參數(shù)確定的，但我們在最后估計(jì)了均值和方差，這個(gè)假設(shè)就被違反了，故而這個(gè)測試統(tǒng)計(jì)量的P值是含偏的，這個(gè)用法是錯(cuò)誤的。分布尾部最后，我們可以檢查分布的右尾，我們可以使用分位點(diǎn)函數(shù)ppf，其為cdf函數(shù)的逆，來獲得臨界值，或者更直

40、接的，我們可以使用殘存函數(shù)的逆來辦。 crit01, crit05, crit10 = stats.t.ppf(1-0.01, 1-0.05, 1-0.10, 10) print critical values from ppf at 1%, 5% and 10% %8.4f %8.4f %8.4f% (crit01, crit05, crit10)critical values from ppf at 1%, 5% and 10% 2.7638 1.8125 1.3722 print critical values from isf at 1%, 5% and 10% %8.4f %8.4f

41、 %8.4f% tuple(stats.t.isf(0.01,0.05,0.10,10)critical values from isf at 1%, 5% and 10% 2.7638 1.8125 1.3722 freq01 = np.sum(xcrit01) / float(n) * 100 freq05 = np.sum(xcrit05) / float(n) * 100 freq10 = np.sum(xcrit10) / float(n) * 100 print sample %-frequency at 1%, 5% and 10% tail %8.4f %8.4f %8.4f%

42、 (freq01, freq05, freq10)sample %-frequency at 1%, 5% and 10% tail 1.4000 5.8000 10.5000在這三種情況中，我們的樣本有有一個(gè)更重的尾部，即實(shí)際在理論分界值右邊的概率要高于理論值。我們可以通過使用更大的樣本來獲得更好的擬合。在以下情況經(jīng)驗(yàn)頻率已經(jīng)很接近理論概率了，但即使我們重復(fù)這個(gè)過程若干次，波動(dòng)依然會保持在這個(gè)程度。 freq05l = np.sum(stats.t.rvs(10, size=10000) crit05) / 10000.0 * 100 print larger sample %-freque

43、ncy at 5% tail %8.4f% freq05llarger sample %-frequency at 5% tail 4.8000我們也可以比較它與正態(tài)分布的尾部，其有一個(gè)輕的多的尾部： print tail prob. of normal at 1%, 5% and 10% %8.4f %8.4f %8.4f% . tuple(stats.norm.sf(crit01, crit05, crit10)*100)tail prob. of normal at 1%, 5% and 10% 0.2857 3.4957 8.5003卡方檢驗(yàn)可以被用來測試，是否一個(gè)有限的分類觀測值頻率

44、與假定的理論概率分布具有顯著差異。 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) print crit -Inf -2.76376946 -1.81246112 -1.37218364 1.37218364 1.81246112 2.76376946 Inf n_sample = x.size freqcount = np.histogram(x, bins=crit)0 tprob = np.diff(quantiles) nprob = np.d

45、iff(stats.norm.cdf(crit) tch, tpval = stats.chisquare(freqcount, tprob*n_sample) nch, npval = stats.chisquare(freqcount, nprob*n_sample) print chisquare for t: chi2 = %6.3f pvalue = %6.4f % (tch, tpval)chisquare for t: chi2 = 2.300 pvalue = 0.8901 print chisquare for normal: chi2 = %6.3f pvalue = %6

46、.4f % (nch, npval)chisquare for normal: chi2 = 64.605 pvalue = 0.0000我們看到當(dāng)t分布檢驗(yàn)沒被拒絕時(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)。 tdof, tloc, tscale = stats.t.fit(x) nloc, nscale = stats.norm.fit(x) tprob = np.diff(stats.t.cdf

47、(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 chisquare for t: chi2 = %6.3f pvalue = %6.4f % (tch, tpval)chisquare for t: chi2 = 1

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

49、skewtest teststat = %6.3f pvalue = %6.4f % stats.skewtest(x)normal skewtest teststat = 2.785 pvalue = 0.0054 print normal kurtosistest teststat = %6.3f pvalue = %6.4f % stats.kurtosistest(x)normal kurtosistest teststat = 4.757 pvalue = 0.0000將這兩個(gè)檢驗(yàn)組合起來的正態(tài)性檢驗(yàn) print normaltest teststat = %6.3f pvalue

50、= %6.4f % stats.normaltest(x)normaltest teststat = 30.379 pvalue = 0.0000在所有三個(gè)測試中，P值是非常低的，所以我們可以拒絕我們的樣本的峰度與偏度與正態(tài)分布相同的假設(shè)。當(dāng)我們的樣本標(biāo)準(zhǔn)化之后，我們依舊得到相同的結(jié)果。 print normaltest teststat = %6.3f pvalue = %6.4f % . stats.normaltest(x-x.mean()/x.std()normaltest teststat = 30.379 pvalue = 0.0000因?yàn)檎龖B(tài)性被很強(qiáng)的拒絕了，所以我們可以檢查這種

51、檢驗(yàn)方式是否可以有效地作用到其他情況中。 print normaltest teststat = %6.3f pvalue = %6.4f % stats.normaltest(stats.t.rvs(10, size=100)normaltest teststat = 4.698 pvalue = 0.0955 print normaltest teststat = %6.3f pvalue = %6.4f % stats.normaltest(stats.norm.rvs(size=1000)normaltest teststat = 0.613 pvalue = 0.7361我們檢驗(yàn)了小樣本的t分布樣本的觀測值以及一個(gè)大樣本的正態(tài)分布觀測值，在這兩種情況中我們都不能拒絕其來自正態(tài)分布的空假設(shè)。得到這樣的結(jié)果是因?yàn)榍罢呤且驗(yàn)闊o法區(qū)分小樣本下的t分布，后者是因?yàn)樗緛砭蛠碜哉龖B(tài)分布。比較兩個(gè)樣本接下來，我們有兩個(gè)