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模塊的特性。這里我們的意圖是提供給使用者一個關(guān)于這個包的實用性知識。我們推薦reference manual來介紹更多的細(xì)節(jié)。

注意：這個文檔還在發(fā)展中。

隨機變量

有一些通用的概率分布類被封裝在continuous random variables以及 discrete random variables中。有80多個連續(xù)性隨機變量(RVs)以及10余個離散隨機變量已經(jīng)用 這些類建立。同樣，新的程序和分布可以被用戶新建(如果你構(gòu)造了一個，請?zhí)峁┧o我們幫助發(fā)展 這個包)。

所有統(tǒng)計函數(shù)被放在子包scipy.stats中，且有這些函數(shù)的一個幾乎完整的列表可以使用 info(stats)獲得。這個列表里的隨機變量也可以從stats子包的docstring中獲得介紹。

在接下來的討論中，我們著重于連續(xù)性隨機變量(RVs)。幾乎所有離散變量也符合下面的討論， 盡管我們將“離散分布的特殊之處”指出它們的一些區(qū)別。

下面的示例代碼我們假設(shè)scipy.stats包已被下述方式導(dǎo)入。

>>> from scipy import stats

有些例子假設(shè)對象被這樣的方式導(dǎo)入(不用輸完整路徑)了。

>>> from scipy.stats import norm

獲得幫助

所有分布可以使用help函數(shù)得到解釋。為獲得這些信息只需要使用像這樣的簡單調(diào)用：

>>> print norm.__doc__

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

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

bounds of distribution lower: -inf, upper: inf

我們可以通過調(diào)用dir(norm)來獲得關(guān)于這個(正態(tài))分布的所有方法和屬性。應(yī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)省獲得所有的可用分布的信息。

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

通用方法

連續(xù)隨機變量的主要公共方法如下：

rvs:隨機變量(就是從這個分布中抽一些樣本)

pdf：概率密度函數(shù)。

cdf：累計分布函數(shù)

sf：殘存函數(shù)(1-CDF)

ppf：分位點函數(shù)(CDF的逆)

isf：逆殘存函數(shù)(sf的逆)

stats:返回均值，方差，(費舍爾)偏態(tài)，(費舍爾)峰度。

moment:分布的非中心矩。

讓我們使用一個標(biāo)準(zhǔn)正態(tài)(normal)隨機變量(RV)作為例子。

>>> norm.cdf(0)

0.5

為了計算在一個點上的cdf，我們可以傳遞一個列表或一個numpy數(shù)組。

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

相應(yīng)的，像pdf,cdf之類的簡單方法可以用np.vectorize進(jìn)行矢量化.

一些其他的實用通用方法:

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

(0.0, 1.0, 1.0)

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

(array(0.0), array(1.0))

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

>>> norm.ppf(0.5)

0.0

為了產(chǎn)生一個隨機變量列,使用size關(guān)鍵字參數(shù)。

>>> norm.rvs(size=5)

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

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

>>> norm.rvs(5)

7.131624370075814

欲知其意，請看下一部分的內(nèi)容。

偏移(Shifting)與縮放(Scaling)

所有連續(xù)分布可以操縱loc以及scale參數(shù)調(diào)整分布的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))

通常經(jīng)標(biāo)準(zhǔn)化的分布的隨機變量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)=1?exp(?λx)

通過像上面那樣使用scale，可以看到如何得到目標(biāo)期望值。

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

最后，聯(lián)系起我們在前面段落中留下的norm.rvs(5)的問題。事實上，像這樣調(diào)用一個分布， 其第一個參數(shù)，像之前的5，是把loc參數(shù)調(diào)到了5，讓我們看：

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

4.983550784784704

在這里，為解釋最后一段的輸出：norm.rvs(5)產(chǎn)生了一個正態(tài)分布變量，其期望，即loc=5.

我傾向于明確的使用loc,scale作為關(guān)鍵字而非像上面那樣依賴參數(shù)的順序。 因為這看上去有點令人困惑。我們在我們解釋“凍結(jié)RV”的主題之前澄清這一點。

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

雖然一般連續(xù)隨機變量都可以通過賦予loc和scale參數(shù)進(jìn)行偏移和縮放，但一些分布還需要 額外的形態(tài)參數(shù)確定其形態(tài)。作為例子，看到這個伽馬分布，這是它的密度函數(shù)

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

它要求一個形態(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.numargs

1

>>> 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)該得到的。

廣播

像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)， 則進(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.68099799])

離散分布的特殊之處

離散分布的方法的大多數(shù)與連續(xù)分布很類似。當(dāng)然像pdf被更換為密度函數(shù)pmf，沒有估計方法， 像fit就不能用了。而scale不是一個合法的關(guān)鍵字參數(shù)。Location參數(shù)， 關(guān)鍵字loc則仍然可以使用用于位移。

cdf的計算要求一些額外的關(guān)注。在連續(xù)分布的情況下，累積分布函數(shù)在大多數(shù)標(biāo)準(zhǔn)情況下是嚴(yán)格遞增的， 所以有唯一的逆。而cdf在離散分布卻一般是階躍函數(shù)，所以cdf的逆，分位點函數(shù)，要求一個不同的定義：

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

為了更多信息可以看這里。

我們可以看這個超幾何分布的例子

>>> from scipy.stats import hypergeom

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

如果我們在一些整數(shù)點使用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的函數(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ù)的極大似然估計，包括location與scale

fit_loc_scale: 給定形態(tài)參數(shù)確定下的location和scale參數(shù)的估計

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

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

性能問題與注意事項

分布方法的性能與運行速度根據(jù)分布的不同表現(xiàn)差異極大。方法的結(jié)果可以由兩種方式獲得， 精確計算或使用獨立于各具體分布的通用算法。

精確計算一般更快。為了進(jìn)行精確計算，要么直接使用解析公式，要么使用scipy.special中的 函數(shù)，對于rvs還可以使用numpy.random里的函數(shù)。

另一方面，如果不能進(jìn)行精確計算，將使用通用方法進(jìn)行計算。于是為了定義一個分布， 只有pdf異或cdf是必須的；通用方法使用數(shù)值積分和求根法進(jìn)行求解。作為例子， rgh = stats.gausshyper.rvs(0.5, 2, 2, 2, size=100)以這種方式創(chuàng)建了100個隨機變量 (抽了100個值)，這在我的電腦上花了19秒(譯者：我花了3.5秒)， 對比取一百萬個標(biāo)準(zhǔn)正態(tài)分布的值只需要1秒。

遺留問題

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

分布在很多參數(shù)區(qū)間上的值被測試過了，然而在一些奇葩的臨界條件，仍然可能有錯誤的值存在。

fit的極大似然估計以默認(rèn)值作為初始參數(shù)將不會工作的很好，用戶必須指派合適的初始參數(shù)。 并且，對于一些分布使用極大似然估計本身就不是一個好的選擇。

構(gòu)造具體的分布

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

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

創(chuàng)建連續(xù)分布是非常簡單的.

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

注意這種用法的性能問題，參見“性能問題與注意事項”一節(jié)。這種缺乏信息的通用計算可能非常慢。 而且再看看下面這個準(zhǔn)確性的例子：

>>> from scipy.integrate import quad

>>> quad(deterministic.pdf, -1e-1, 1e-1)

(4.163336342344337e-13, 0.0)

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

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

(1.000076872229173, 0.0010625571718182458)

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

繼承rv_discrete類

在之后我們使用stats.rv_discrete產(chǎn)生一個離散分布，其有一個整數(shù)區(qū)間截斷概率。

通用信息

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

>>> from scipy.stats import rv_discrete

>>> help(rv_discrete)

從中我們得知：

“你可以構(gòu)建任意一個像P(X=xk)=pk一樣形式的離散rv，通過傳遞(xk,pk)元組序列給 rv_discrete初始化方法(通過values=keyword方式)，但其不能有0概率值。”

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

keyword必須給出。

Xk必須是整數(shù)

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

事實上，如果最后兩個要求沒有被滿足，一個異常將被拋出或者導(dǎo)致一個錯誤的數(shù)值。

一個例子

讓我們開始辦，首先

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

現(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) # 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分布。這也是在驗證 是否隨機數(shù)是以正確的方式產(chǎn)生的。

卡方測試要求起碼在每個子區(qū)間(bin)里具有最小數(shù)目的觀測值。我們組合末端子區(qū)間進(jìn)大子區(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 'chisquare for normdiscrete: chi2 = %6.3f pvalue = %6.4f' % (ch2, pval)

chisquare for normdiscrete: chi2 = 12.466 pvalue = 0.4090

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

樣本分析

首先，我們創(chuàng)建一些隨機變量。我們設(shè)置一個種子所以每次我們都可以得到相同的結(jié)果以便觀察。 作為一個例子，我們從t分布中抽一個樣本。

>>> np.random.seed(282629734)

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

這里，我們設(shè)置了t分布的形態(tài)參數(shù)，在這里就是自由度，設(shè)為10.使用size=1000表示 我們的樣本由1000個抽樣是獨立的(偽)。當(dāng)我們不指派loc和scale時，它們具有默認(rèn)值0和1.

描述統(tǒng)計

X是一個numpy數(shù)組。我們可以直接調(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)

0.0140610663985 1.28899386208

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

>>> 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檢驗是否樣本與給定均值(這里是理論均值)存在統(tǒng)計顯著差異。

>>> 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”這個假設(shè)， 0是標(biāo)準(zhǔn)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檢驗)被被用來檢驗樣本是否來自一個標(biāo)準(zhǔn)的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分布的假設(shè)。在實際應(yīng)用中， 我們不能知道潛在的分布到底是什么。如果我們使用KS檢驗我們的樣本對照正態(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時。如果我們標(biāo)準(zhǔn)化我們的樣本并且 測試它比照正態(tài)分布，那么P值將又一次很高我們將還是不能拒絕假設(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檢驗假設(shè)我們比照的分布就是以給定的參數(shù)確定的，但我們在最后估計了均值和方差， 這個假設(shè)就被違反了，故而這個測試統(tǒng)計量的P值是含偏的，這個用法是錯誤的。

分布尾部

最后，我們可以檢查分布的右尾，我們可以使用分位點函數(shù)ppf，其為cdf函數(shù)的逆，來獲得臨界值， 或者更直接的，我們可以使用殘存函數(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 %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

在這三種情況中，我們的樣本有有一個更重的尾部，即實際在理論分界值右邊的概率要高于理論值。 我們可以通過使用更大的樣本來獲得更好的擬合。在以下情況經(jīng)驗頻率已經(jīng)很接近理論概率了， 但即使我們重復(fù)這個過程若干次，波動依然會保持在這個程度。

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

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

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

我們看到當(dāng)t分布檢驗沒被拒絕時標(biāo)準(zhǔn)正態(tài)分布卻被完全拒絕。在我們的樣本區(qū)分出這兩個分布后， 我們可以先進(jìn)行擬合確定scale與location再檢查擬合后的分布的差異性。

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

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

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

正態(tài)分布的特殊檢驗

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

首先我們檢驗分布的峰度和偏度是否顯著地與正態(tài)分布的對應(yīng)值相差異。

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

將這兩個檢驗組合起來的正態(tài)性檢驗

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

normaltest teststat = 30.379 pvalue = 0.0000

在所有三個測試中，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

因為正態(tài)性被很強的拒絕了，所以我們可以檢查這種檢驗方式是否可以有效地作用到其他情況中。

>>> 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ài)分布觀測值，在這兩種情況中我們 都不能拒絕其來自正態(tài)分布的空假設(shè)。得到這樣的結(jié)果是因為前者是因為無法區(qū)分小樣本下的t分布， 后者是因為它本來就來自正態(tài)分布。

比較兩個樣本

接下來，我們有兩個分布，其可以判定為相同或者來自不同的分布，以及我們希望測試是否這些 樣本有相同的統(tǒng)計特征。

均值

以相同的均值產(chǎn)生的樣本進(jì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.54890361750888583, 0.5831943748663857)

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

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

>>> stats.ttest_ind(rvs1, rvs3)

(-4.5334142901750321, 6.507128186505895e-006)

對于兩個不同的樣本進(jìn)行的KS檢驗

在這個例子中我們使用兩個同分布的樣本進(jìn)行檢驗.設(shè)因為P值很高,毫不奇怪我們不能拒絕原假設(shè)。

>>> stats.ks_2samp(rvs1, rvs2)

(0.025999999999999995, 0.99541195173064878)

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

>>> stats.ks_2samp(rvs1, rvs3)

(0.11399999999999999, 0.0027132103661283141)

核密度估計

一個常見的統(tǒng)計學(xué)問題是從一個樣本中估計隨機變量的概率密度分布函數(shù)(PDF) 這個問題被稱為密度估計，對此最著名的工具是直方圖。直方圖是一個很好的可視化工具 (主要是因為每個人都理解它)。但是對于對于數(shù)據(jù)特征的利用卻并不是非常有效率。

核密度估計(KDE對于這個問題)是一個更有效的工具。這個gaussian_kde估計方法可以被用來估計 單元或多元數(shù)據(jù)的PDF。它在數(shù)據(jù)呈單峰的時候工作的最好，但也可以在多峰情況下工作。

單元估計

我們以一個最小數(shù)據(jù)集來觀察gaussian_kde是如何工作的，以及帶寬(bandwidth)的不同選擇方式。 PDF對應(yīng)的數(shù)據(jù)被以藍(lán)線的形式顯示在圖像的底端(被稱為毯圖(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規(guī)則以及Silverman規(guī)則下的結(jié)果幾乎沒有差異。以及帶寬的選擇相比較于 數(shù)據(jù)的稀少顯得太寬。我們可以定義我們的帶寬函數(shù)以獲得一個更少平滑的結(jié)果。

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

我們看到如果我們設(shè)置帶寬為非常狹窄，則獲得PDF的估計退化為圍繞在數(shù)據(jù)點的簡單的高斯和。

我們現(xiàn)在使用更真實的例子，并且看看在兩種帶寬選擇規(guī)則中的差異。這些規(guī)則被認(rèn)為在 正態(tài)分布上很好用，但即使是偏離正態(tài)的單峰分布上它也工作的很好。作為一個非正態(tài)分布， 我們采用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()

下面我們看到這個一個寬一個窄的雙峰分布。可以想到結(jié)果將難達(dá)到以十分近似， 因為每個峰需要不同的帶寬去擬合。

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

正如預(yù)想的，KDE并沒有很好的趨近正確的PDF，因為雙峰分布的形狀不同。通過使用默認(rèn)帶寬 (Scott*0.5)我們可以做得更好，再使用更小的帶寬將使平滑性受到影響。這里我們真正需要 的是非均勻(自適應(yīng))帶寬。

多元估計

通過gaussian_kde我們可以像單元估計那樣進(jìn)行多元估計。我們現(xiàn)在來解決二元情況， 首先我們產(chǎn)生一些隨機的二元數(shù)據(jù)。

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

然后我們對這些數(shù)據(jù)使用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形式顯示出來，并且在上面畫出每個數(shù)據(jù)點。

>>> 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位技術(shù)專家面對面20年技術(shù)見證，附贈技術(shù)全景圖

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。