scipy模塊stats文檔

https://github.com/yiyuezhuo/scipy.stats-doc-ch

https://docs.scipy.org/doc/scipy/reference/tutorial/stats.html

http://blog.csdn.net/pipisorry/article/details/49515215

 

介紹

在這個教程咱們討論一部分scipy.stats模塊的特性。這裏咱們的意圖是提供給使用者一個關於這個包的實用性知識。咱們推薦reference manual來介紹更多的細節。

注意：這個文檔還在發展中。

 

隨機變量

有一些通用的機率分佈類被封裝在continuous random variables以及 discrete random variables中。有80多個連續性隨機變量(RVs)以及10餘個離散隨機變量已經用 這些類創建。一樣，新的程序和分佈能夠被用戶新建（若是你構造了一個，請提供它給咱們幫助發展 這個包）。

全部統計函數被放在子包scipy.stats中，且有這些函數的一個幾乎完整的列表可使用 info(stats)得到。這個列表裏的隨機變量也能夠從stats子包的docstring中得到介紹。

在接下來的討論中，咱們着重於連續性隨機變量(RVs)。幾乎全部離散變量也符合下面的討論， 儘管咱們將「離散分佈的特殊之處」指出它們的一些區別。

下面的示例代碼咱們假設scipy.stats包已被下述方式導入。

>>> from scipy import stats

有些例子假設對象被這樣的方式導入（不用輸完整路徑）了。

>>> from scipy.stats import norm

得到幫助

全部分佈可使用help函數獲得解釋。爲得到這些信息只須要使用像這樣的簡單調用：

>>> print norm.__doc__

做爲例子，咱們用這種方式獲取分佈的上下界

>>> print 'bounds of distribution lower: %s, upper: %s' % (norm.a,norm.b) bounds of distribution lower: -inf, upper: inf

咱們能夠經過調用dir(norm)來得到關於這個（正態）分佈的全部方法和屬性。應該看到， 一些方法是私有方法儘管其並無以名稱表示出來（好比它們前面沒有如下劃線開頭）， 好比veccdf就只用於內部計算（試圖使用那些方法將引起警告，由於它們可能會在後續開發中被移除）

爲了得到真正的主要方法，咱們列舉凍結分佈的方法（咱們將在下文解釋何謂凍結）

>>> 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']

最後，咱們能經過內省得到全部的可用分佈的信息。

>>> import warnings >>> warnings.simplefilter('ignore', DeprecationWarning) >>> dist_continu = [d for d in dir(stats) if ... isinstance(getattr(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 discrete distributions: 12

通用方法

連續隨機變量的主要公共方法以下：

• rvs:隨機變量（就是從這個分佈中抽一些樣本）
• pdf：機率密度函數。
• cdf：累計分佈函數
• sf：殘存函數（1-CDF）
• ppf：分位點函數（CDF的逆）
• isf：逆殘存函數（sf的逆）
• stats:返回均值，方差，（費舍爾）偏態，（費舍爾）峯度。
• moment:分佈的非中心矩。

讓咱們使用一個標準正態(normal)隨機變量(RV)做爲例子。

>>> norm.cdf(0) 0.5

爲了計算在一個點上的cdf，咱們能夠傳遞一個列表或一個numpy數組。

>>> norm.cdf([-1., 0, 1]) array([ 0.15865525, 0.5 , 0.84134475]) >>> import numpy as np >>> norm.cdf(np.array([-1., 0, 1])) array([ 0.15865525, 0.5 , 0.84134475])

相應的，像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))

爲了找到一個分佈的中中位數，咱們可使用分位數函數ppf，它是cdf的逆。

>>> norm.ppf(0.5) 0.0

爲了產生一個隨機變量列,使用size關鍵字參數。

>>> norm.rvs(size=5) array([-0.35687759, 1.34347647, -0.11710531, -1.00725181, -0.51275702])

不要認爲norm.rvs(5)產生了五個變量。

>>> norm.rvs(5) 7.131624370075814

欲知其意，請看下一部分的內容。

偏移(Shifting)與縮放(Scaling)

全部連續分佈能夠操縱loc以及scale參數調整分佈的location和scale屬性。做爲例子， 標準正態分佈的location是均值而scale是標準差。

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

一般經標準化的分佈的隨機變量X能夠經過變換(X-loc)/scale得到。它們的默認值是loc=0以及scale=1.

聰明的使用loc與scale能夠幫助以靈活的方式調整標準分佈達到所想要的效果。 爲了進一步說明縮放的效果，下面給出指望爲1/λ指數分佈的cdf。

F(x)=1−exp(−λ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 = 1, scale = 4) array([ 0. , 0. , 0.25, 0.5 , 0.75, 1. ])

最後，聯繫起咱們在前面段落中留下的norm.rvs(5)的問題。事實上，像這樣調用一個分佈， 其第一個參數，像以前的5，是把loc參數調到了5，讓咱們看：

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

在這裏，爲解釋最後一段的輸出：norm.rvs(5)產生了一個正態分佈變量，其指望，即loc=5.

我傾向於明確的使用loc,scale做爲關鍵字而非像上面那樣依賴參數的順序。 由於這看上去有點使人困惑。咱們在咱們解釋「凍結RV」的主題以前澄清這一點。

形態(shape)參數

雖然通常連續隨機變量均可以經過賦予loc和scale參數進行偏移和縮放，但一些分佈還須要 額外的形態參數肯定其形態。做爲例子，看到這個伽馬分佈，這是它的密度函數

γ(x,a)=λ(λx)a−1Γ(a)e−λx,

它要求一個形態參數a。注意到λ的設置能夠經過設置scale關鍵字爲1/λ進行。

讓咱們檢查伽馬分佈的形態參數的名字的數量。（咱們從上面知道其應該爲1）

>>> from scipy.stats import gamma >>> gamma.numargs 1 >>> gamma.shapes 'a'

如今咱們設置形態變量的值爲1以變成指數分佈。因此咱們能夠容易的比較是否獲得了咱們所指望的結果。

>>> gamma(1, scale=2.).stats(moments="mv") (array(2.0), array(4.0))

注意咱們也能夠以關鍵字的方式指定形態參數：

>>> gamma(a=1, scale=2.).stats(moments="mv") (array(2.0), array(4.0))

凍結分佈

不斷地傳遞loc與scale關鍵字最終會讓人厭煩。而凍結RV的概念被用來解決這個問題。

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

經過使用rv，在任何狀況下咱們再也不須要包含scale與形態參數。顯然，分佈能夠被多種方式使用， 咱們能夠經過傳遞全部分佈參數給對方法的每次調用（像咱們以前作的那樣）或者能夠對一個分 布對象先凍結參數。讓咱們看看是怎麼回事：

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

這正是咱們應該獲得的。

廣播

像pdf這樣的簡單方法知足numpy的廣播規則。做爲例子，咱們能夠計算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爲自由度的臨界值。因此， 廣播規則與下面調用了兩次isf產生的結果相同。

>>> 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])

可是若是機率數組，如[0.1,0.05,0.01]與自由度數組,如[10,11,12]具備相同的數組形態， 則進行對應匹配，咱們能夠分別獲得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])

離散分佈的特殊之處

離散分佈的方法的大多數與連續分佈很相似。固然像pdf被更換爲密度函數pmf，沒有估計方法， 像fit就不能用了。而scale不是一個合法的關鍵字參數。Location參數， 關鍵字loc則仍然可使用用於位移。

cdf的計算要求一些額外的關注。在連續分佈的狀況下，累積分佈函數在大多數標準狀況下是嚴格遞增的， 因此有惟一的逆。而cdf在離散分佈卻通常是階躍函數，因此cdf的逆，分位點函數，要求一個不一樣的定義：

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

爲了更多信息能夠看這裏。

咱們能夠看這個超幾何分佈的例子

>>> from scipy.stats import hypergeom >>> [M, n, N] = [20, 7, 12]

若是咱們在一些整數點使用cdf，則它們的cdf值再做用ppf會回到開始的值。

>>> x = np.arange(4)*2 >>> x array([0, 2, 4, 6]) >>> prb = hypergeom.cdf(x, M, n, N) >>> prb array([ 0.0001031991744066, 0.0521155830753351, 0.6083591331269301, 0.9897832817337386]) >>> hypergeom.ppf(prb, M, n, N) array([ 0., 2., 4., 6.])

若是咱們使用的值不是cdf的函數值，則咱們獲得一個更高的值。

>>> 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.])

分佈擬合

非凍結分佈的參數估計的主要方法：

• fit：分佈參數的極大似然估計，包括location與scale
• fit_loc_scale: 給定形態參數肯定下的location和scale參數的估計
• nnlf:負對數似然函數
• expect:計算函數pdf或pmf的指望值。

性能問題與注意事項

分佈方法的性能與運行速度根據分佈的不一樣表現差別極大。方法的結果能夠由兩種方式得到， 精確計算或使用獨立於各具體分佈的通用算法。

精確計算通常更快。爲了進行精確計算，要麼直接使用解析公式，要麼使用scipy.special中的 函數，對於rvs還可使用numpy.random裏的函數。

另外一方面，若是不能進行精確計算，將使用通用方法進行計算。因而爲了定義一個分佈， 只有pdf異或cdf是必須的；通用方法使用數值積分和求根法進行求解。做爲例子， rgh = stats.gausshyper.rvs(0.5, 2, 2, 2, size=100)以這種方式建立了100個隨機變量 （抽了100個值），這在個人電腦上花了19秒（譯者：我花了3.5秒）， 對比取一百萬個標準正態分佈的值只須要1秒。

遺留問題

scipy.stats裏的分佈最近進行了升級而且被仔細的檢查過了，不過仍有一些問題存在。

• 分佈在不少參數區間上的值被測試過了，然而在一些奇葩的臨界條件，仍然可能有錯誤的值存在。
• fit的極大似然估計以默認值做爲初始參數將不會工做的很好，用戶必須指派合適的初始參數。 而且，對於一些分佈使用極大似然估計自己就不是一個好的選擇。

構造具體的分佈

下一個例子展現瞭如何創建你本身的分佈。更多的例子見分佈用法以及統計檢驗

建立一個連續分佈，繼承rv_continuous類

建立連續分佈是很是簡單的.

>>> from scipy import stats >>> class deterministic_gen(stats.rv_continuous): ... def _cdf(self, x): ... return np.where(x < 0, 0., 1.) ... def _stats(self): ... return 0., 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])

注意這種用法的性能問題，參見「性能問題與注意事項」一節。這種缺少信息的通用計算可能很是慢。 並且再看看下面這個準確性的例子：

>>> from scipy.integrate import quad >>> quad(deterministic.pdf, -1e-1, 1e-1) (4.163336342344337e-13, 0.0)

但這並非對pdf積分的正確的結果，實際上結果應爲1.讓咱們將積分變得更小一些。

>>> quad(deterministic.pdf, -1e-3, 1e-3) # warning removed (1.000076872229173, 0.0010625571718182458)

這樣看上去好多了，然而，問題自己來源於pdf不是來自包給定的類的定義。

繼承rv_discrete類

在以後咱們使用stats.rv_discrete產生一個離散分佈，其有一個整數區間截斷機率。

通用信息

通用信息能夠從 rv_discrete的 docstring中獲得。

>>> from scipy.stats import rv_discrete >>> help(rv_discrete)

從中咱們得知：

「你能夠構建任意一個像P(X=xk)=pk同樣形式的離散rv，經過傳遞(xk,pk)元組序列給 rv_discrete初始化方法(經過values=keyword方式)，但其不能有0機率值。」

接下來，還有一些進一步的要求：

• keyword必須給出。
• Xk必須是整數
• 小數的有效位數應當被給出。

事實上，若是最後兩個要求沒有被知足，一個異常將被拋出或者致使一個錯誤的數值。

一個例子

讓咱們開始辦，首先

>>> npoints = 20 # number of integer support points of the distribution minus 1 >>> npointsh = npoints / 2 >>> npointsf = float(npoints) >>> nbound = 4 # 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 >>> 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')

如今咱們已經定義了這個分佈，咱們能夠調用其全部常規的離散分佈方法。

>>> 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'))

測試上面的結果

讓咱們產生一個隨機樣本而且比較連續機率的狀況。

>>> n_sample = 500 >>> np.random.seed(87655678) # fix the seed for replicability >>> rvs = normdiscrete.rvs(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 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.00000000e+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.62125309e+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]]

接下來，咱們能夠測試，是否咱們的樣本取自於一個normdiscrete分佈。這也是在驗證 是否隨機數是以正確的方式產生的。

卡方測試要求起碼在每一個子區間（bin）裏具備最小數目的觀測值。咱們組合末端子區間進大子區間 因此它們如今包含了足夠數量的觀測值。

>>> 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 'chisquare for normdiscrete: chi2 = %6.3f pvalue = %6.4f' % (ch2, pval) chisquare for normdiscrete: chi2 = 12.466 pvalue = 0.4090

P值在這個狀況下是不顯著地，因此咱們能夠斷言咱們的隨機樣本的確是由此分佈產生的。

樣本分析

首先，咱們建立一些隨機變量。咱們設置一個種子因此每次咱們均可以獲得相同的結果以便觀察。 做爲一個例子，咱們從t分佈中抽一個樣本。

>>> np.random.seed(282629734) >>> x = stats.t.rvs(10, size=1000)

這裏，咱們設置了t分佈的形態參數，在這裏就是自由度，設爲10.使用size=1000表示 咱們的樣本由1000個抽樣是獨立的（僞）。當咱們不指派loc和scale時，它們具備默認值0和1.

描述統計

X是一個numpy數組。咱們能夠直接調用它的方法。

>>> 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) 0.0140610663985 1.28899386208

如何比較分佈自己和它的樣本的指標？

>>> 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' >>> 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.describe用的是無偏的方差估計量，而np.var卻用的是有偏的估計量。

T檢驗和KS檢驗

咱們可使用t檢驗是否樣本與給定均值（這裏是理論均值）存在統計顯著差別。

>>> print 't-statistic = %6.3f pvalue = %6.4f' % stats.ttest_1samp(x, m) t-statistic = 0.391 pvalue = 0.6955

P值爲0.7，這表明第一類錯誤的機率，在例子中，爲10%。咱們不能拒絕「該樣本均值爲0」這個假設， 0是標準t分佈的理論均值。

>>> tt = (sm-m)/np.sqrt(sv/float(n)) # t-statistic for mean >>> pval = 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檢驗（KS檢驗）被被用來檢驗樣本是否來自一個標準的t分佈。

>>> print 'KS-statistic D = %6.3f pvalue = %6.4f' % stats.kstest(x, 't', (10,)) KS-statistic D = 0.016 pvalue = 0.9606

又一次獲得了很高的P值。因此咱們不能拒絕樣本是來自t分佈的假設。在實際應用中， 咱們不能知道潛在的分佈究竟是什麼。若是咱們使用KS檢驗咱們的樣本對照正態分佈， 那麼咱們將也不能拒絕咱們的樣本是來自正態分佈的，在這種狀況下P值爲0.4左右。

>>> print 'KS-statistic D = %6.3f pvalue = %6.4f' % stats.kstest(x,'norm') KS-statistic D = 0.028 pvalue = 0.3949

不管如何，標準正態分佈有1的方差，當咱們的樣本有1.29時。若是咱們標準化咱們的樣本而且 測試它比照正態分佈，那麼P值將又一次很高咱們將仍是不能拒絕假設是來自正態分佈的。

>>> 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檢驗假設咱們比照的分佈就是以給定的參數肯定的，但咱們在最後估計了均值和方差， 這個假設就被違反了，故而這個測試統計量的P值是含偏的，這個用法是錯誤的。

分佈尾部

最後，咱們能夠檢查分佈的右尾，咱們可使用分位點函數ppf，其爲cdf函數的逆，來得到臨界值， 或者更直接的，咱們可使用殘存函數的逆來辦。

>>> 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 %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(x>crit01) / float(n) * 100 >>> freq05 = np.sum(x>crit05) / float(n) * 100 >>> freq10 = np.sum(x>crit10) / float(n) * 100 >>> print 'sample %%-frequency at 1%%, 5%% and 10%% tail %8.4f %8.4f %8.4f'% (freq01, freq05, freq10) sample %-frequency at 1%, 5% and 10% tail 1.4000 5.8000 10.5000

在這三種狀況中，咱們的樣本有有一個更重的尾部，即實際在理論分界值右邊的機率要高於理論值。 咱們能夠經過使用更大的樣原本得到更好的擬合。在如下狀況經驗頻率已經很接近理論機率了， 但即便咱們重複這個過程若干次，波動依然會保持在這個程度。

>>> freq05l = np.sum(stats.t.rvs(10, size=10000) > crit05) / 10000.0 * 100 >>> print 'larger sample %%-frequency at 5%% tail %8.4f'% freq05l larger sample %-frequency at 5% tail 4.8000

咱們也能夠比較它與正態分佈的尾部，其有一個輕的多的尾部：

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

卡方檢驗能夠被用來測試，是否一個有限的分類觀測值頻率與假定的理論機率分佈具備顯著差別。

>>> 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.diff(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.4f' % (nch, npval) chisquare for normal: chi2 = 64.605 pvalue = 0.0000

咱們看到當t分佈檢驗沒被拒絕時標準正態分佈卻被徹底拒絕。在咱們的樣本區分出這兩個分佈後， 咱們能夠先進行擬合肯定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 'chisquare for t: chi2 = %6.3f pvalue = %6.4f' % (tch, tpval) chisquare for t: chi2 = 1.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

在通過參數調整以後，咱們仍然能夠以5%水平拒絕正態分佈假設。然而卻以95%的p值顯然的不能拒絕t分佈。

正態分佈的特殊檢驗

自從正態分佈變爲統計學中最多見的分佈，就出現了大量的方法用來檢驗一個樣本 是否能夠被當作是來自正態分佈的。

首先咱們檢驗分佈的峯度和偏度是否顯著地與正態分佈的對應值相差別。

>>> print 'normal 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

將這兩個檢驗組合起來的正態性檢驗

>>> print 'normaltest teststat = %6.3f pvalue = %6.4f' % stats.normaltest(x) normaltest teststat = 30.379 pvalue = 0.0000

在全部三個測試中，P值是很是低的，因此咱們能夠拒絕咱們的樣本的峯度與偏度與正態分佈相同的假設。

當咱們的樣本標準化以後，咱們依舊獲得相同的結果。

>>> print 'normaltest teststat = %6.3f pvalue = %6.4f' % \ ... stats.normaltest((x-x.mean())/x.std()) normaltest teststat = 30.379 pvalue = 0.0000

由於正態性被很強的拒絕了，因此咱們能夠檢查這種檢驗方式是否能夠有效地做用到其餘狀況中。

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

咱們檢驗了小樣本的t分佈樣本的觀測值以及一個大樣本的正態分佈觀測值，在這兩種狀況中咱們 都不能拒絕其來自正態分佈的空假設。獲得這樣的結果是由於前者是由於沒法區分小樣本下的t分佈， 後者是由於它原本就來自正態分佈。

比較兩個樣本

接下來，咱們有兩個分佈，其能夠斷定爲相同或者來自不一樣的分佈，以及咱們但願測試是否這些 樣本有相同的統計特徵。

均值

以相同的均值產生的樣本進行檢驗：

>>> 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.54890361750888583, 0.5831943748663857)

以不一樣的均值產生的樣本進行檢驗：

>>> rvs3 = stats.norm.rvs(loc=8, scale=10, size=500) >>> stats.ttest_ind(rvs1, rvs3) (-4.5334142901750321, 6.507128186505895e-006)

對於兩個不一樣的樣本進行的KS檢驗

在這個例子中咱們使用兩個同分布的樣本進行檢驗.設由於P值很高,絕不奇怪咱們不能拒絕原假設。

>>> stats.ks_2samp(rvs1, rvs2) (0.025999999999999995, 0.99541195173064878)

在第二個例子中，因爲均值不一樣，因此咱們能夠拒絕空假設，由P值小於1%。

>>> stats.ks_2samp(rvs1, rvs3) (0.11399999999999999, 0.0027132103661283141)

核密度估計

一個常見的統計學問題是從一個樣本中估計隨機變量的機率密度分佈函數（PDF） 這個問題被稱爲密度估計，對此最著名的工具是直方圖。直方圖是一個很好的可視化工具 （主要是由於每一個人都理解它）。可是對於對於數據特徵的利用卻並非很是有效率。

核密度估計(KDE對於這個問題)是一個更有效的工具。這個gaussian_kde估計方法能夠被用來估計 單元或多元數據的PDF。它在數據呈單峯的時候工做的最好，但也能夠在多峯狀況下工做。

單元估計

咱們以一個最小數據集來觀察gaussian_kde是如何工做的，以及帶寬（bandwidth）的不一樣選擇方式。 PDF對應的數據被以藍線的形式顯示在圖像的底端（被稱爲毯圖（rug plot））

>>> from scipy import stats >>> import matplotlib.pyplot as plt >>> 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) # rug plot >>> x_eval = np.linspace(-10, 10, num=200) >>> ax.plot(x_eval, kde1(x_eval), 'k-', label="Scott's Rule") >>> ax.plot(x_eval, kde1(x_eval), 'r-', label="Silverman's Rule") >>> plt.show()

咱們看到在Scott規則以及Silverman規則下的結果幾乎沒有差別。以及帶寬的選擇相比較於 數據的稀少顯得太寬。咱們能夠定義咱們的帶寬函數以得到一個更少平滑的結果。

>>> def my_kde_bandwidth(obj, fac=1./5): ... """We use Scott's Rule, multiplied by a constant factor.""" ... return np.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) # rug plot >>> kde3 = stats.gaussian_kde(x1, bw_method=my_kde_bandwidth) >>> ax.plot(x_eval, kde3(x_eval), 'g-', label="With smaller BW") >>> plt.show()

咱們看到若是咱們設置帶寬爲很是狹窄，則得到PDF的估計退化爲圍繞在數據點的簡單的高斯和。

咱們如今使用更真實的例子，而且看看在兩種帶寬選擇規則中的差別。這些規則被認爲在 正態分佈上很好用，但即便是偏離正態的單峯分佈上它也工做的很好。做爲一個非正態分佈， 咱們採用5自由度的t分佈。

import numpy as np import matplotlib.pyplot as plt from scipy import stats np.random.seed(12456) x1 = np.random.normal(size=200) # random data, normal distribution xs = 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) # rug plot ax1.plot(xs, kde1(xs), 'k-', label="Scott's Rule") ax1.plot(xs, kde2(xs), 'b-', label="Silverman's Rule") ax1.plot(xs, stats.norm.pdf(xs), 'r--', label="True PDF") ax1.set_xlabel('x') ax1.set_ylabel('Density') ax1.set_title("Normal (top) and Student's T$_{df=5}$ (bottom) distributions") ax1.legend(loc=1) x2 = stats.t.rvs(5, size=200) # random data, T distribution 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) # rug plot ax2.plot(xs, kde3(xs), 'k-', label="Scott's Rule") ax2.plot(xs, kde4(xs), 'b-', label="Silverman's Rule") ax2.plot(xs, stats.t.pdf(xs, 5), 'r--', label="True PDF") ax2.set_xlabel('x') ax2.set_ylabel('Density') plt.show()

下面咱們看到這個一個寬一個窄的雙峯分佈。能夠想到結果將難達到以十分近似， 由於每一個峯須要不一樣的帶寬去擬合。

>>> from functools import partial >>> loc1, scale1, size1 = (-2, 1, 175) >>> loc2, scale2, size2 = (2, 0.2, 50) >>> x2 = np.concatenate([np.random.normal(loc=loc1, scale=scale1, size=size1), ... np.random.normal(loc=loc2, scale=scale2, size=size2)]) >>> x_eval = np.linspace(x2.min() - 1, x2.max() + 1, 500) >>> kde = stats.gaussian_kde(x2) >>> kde2 = stats.gaussian_kde(x2, bw_method='silverman') >>> kde3 = stats.gaussian_kde(x2, bw_method=partial(my_kde_bandwidth, fac=0.2)) >>> kde4 = stats.gaussian_kde(x2, bw_method=partial(my_kde_bandwidth, fac=0.5)) >>> pdf = stats.norm.pdf >>> bimodal_pdf = pdf(x_eval, loc=loc1, scale=scale1) * float(size1) / x2.size + \ ... pdf(x_eval, loc=loc2, scale=scale2) * float(size2) / x2.size >>> fig = plt.figure(figsize=(8, 6)) >>> ax = fig.add_subplot(111) >>> ax.plot(x2, np.zeros(x2.shape), 'b+', ms=12) >>> ax.plot(x_eval, kde(x_eval), 'k-', label="Scott's Rule") >>> ax.plot(x_eval, kde2(x_eval), 'b-', label="Silverman's Rule") >>> ax.plot(x_eval, kde3(x_eval), 'g-', label="Scott * 0.2") >>> ax.plot(x_eval, kde4(x_eval), 'c-', label="Scott * 0.5") >>> ax.plot(x_eval, bimodal_pdf, 'r--', label="Actual PDF") >>> ax.set_xlim([x_eval.min(), x_eval.max()]) >>> ax.legend(loc=2) >>> ax.set_xlabel('x') >>> ax.set_ylabel('Density') >>> plt.show()

正如預想的，KDE並無很好的趨近正確的PDF，由於雙峯分佈的形狀不一樣。經過使用默認帶寬 （Scott*0.5）咱們能夠作得更好，再使用更小的帶寬將使平滑性受到影響。這裏咱們真正須要 的是非均勻（自適應）帶寬。

多元估計

經過gaussian_kde咱們能夠像單元估計那樣進行多元估計。咱們如今來解決二元狀況， 首先咱們產生一些隨機的二元數據。

>>> def measure(n): ... """Measurement model, return two coupled measurements.""" ... m1 = np.random.normal(size=n) ... m2 = np.random.normal(scale=0.5, size=n) ... return m1+m2, m1-m2 >>> m1, m2 = measure(2000) >>> xmin = m1.min() >>> xmax = m1.max() >>> ymin = m2.min() >>> ymax = m2.max()

而後咱們對這些數據使用KDE：

>>> X, Y = np.mgrid[xmin:xmax:100j, ymin:ymax:100j] >>> positions = np.vstack([X.ravel(), Y.ravel()]) >>> values = np.vstack([m1, m2]) >>> kernel = stats.gaussian_kde(values) >>> Z = np.reshape(kernel.evaluate(positions).T, X.shape)

最後咱們把估計的雙峯分佈以colormap形式顯示出來，而且在上面畫出每一個數據點。

>>> fig = plt.figure(figsize=(8, 6)) >>> ax = fig.add_subplot(111) >>> ax.imshow(np.rot90(Z), cmap=plt.cm.gist_earth_r, ... extent=[xmin, xmax, ymin, ymax]) >>> ax.plot(m1, m2, 'k.', markersize=2) >>> ax.set_xlim([xmin, xmax]) >>> ax.set_ylim([ymin, ymax]) >>> plt.show()