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

與50位技術專家面對面20年技術見證，附贈技術全景圖

總結

以上是生活随笔為你收集整理的python scipy模块文档_scipy模块stats文档的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。