gamma分布密度函數(shù)介紹（Gamma函數(shù)可視化什么意思）-環(huán)球速看

1. Gamma 分布

(資料圖)

1.1 首先要了解一下Gamma 函數(shù)。

Gamma 函數(shù)在實數(shù)域可以表示為：

Γ ( x ) = ∫ 0 ∞ t x − 1 e − t d t \Gamma(x)=\int_0^{\infty} t^{x-1}e^{-t}dtΓ(x)=∫0∞?tx−1e−tdt;

在整數(shù)域可以表示為：

Γ ( n ) = ( n − 1 ) ! \Gamma(n)=(n-1)!Γ(n)=(n−1)!。

此外，Gamma 函數(shù)具有如下性質(zhì)：

Γ ( x + 1 ) = x Γ ( x ) \Gamma(x+1)=x\Gamma(x)Γ(x+1)=xΓ(x)

Γ ( x + 1 ) = x ! \Gamma(x+1)=x!Γ(x+1)=x!

Γ ( 1 ) = 1 \Gamma(1)=1Γ(1)=1

Γ ( 1 2 ) = π \Gamma(\frac{1}{2})=\sqrt{\pi}Γ(21?)=πundefined?

Γ \GammaΓ函數(shù)是階乘在實數(shù)上的推廣。

1.2 Gamma函數(shù)可視化

"""

Beta 分布

@date 2019-05-09

@author wangxianpeng

"""

import numpy as np

from scipy.special import gamma

import matplotlib.pyplot as plt

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

# Gamma函數(shù)(Gamma function)

x = np.linspace(-1, 10, 1000)

plt.plot(x, gamma(x), ls='-', c='k', label='$\Gamma(x)$')

# (x-1)! for x = 1, 2, ..., 6

x2 = np.linspace(1,6,6)

y2 = np.array([1, 1, 2, 6, 24, 120])

plt.plot(x2, y2, marker='*', markersize=12, markeredgecolor='r',

markerfacecolor='r', ls='',c='r', label='$(x-1)!$')

plt.title('Gamma Function')

plt.ylim(0,25)

plt.xlim(0, 6)

plt.xlabel('$x$')

plt.ylabel('$\Gamma(x)$')

plt.legend()

plt.show()

結(jié)果顯示如下：

當(dāng)然這里只取了大于0的部分。

1.3 Gamma 分布

對Gamma函數(shù)做個變形，可以得到如下式子：∫ 0 ∞ t α − 1 e − t d t Γ ( α ) = 1 ? ( 1.3.1 ) \int_0^{\infty} \frac{t^{\alpha-1}e^{-t}dt}{\Gamma(\alpha)}=1 \: (1.3.1)∫0∞?Γ(α)tα−1e−tdt?=1(1.3.1)

取等式左邊積分中的函數(shù)作為概率密度，就得到一個簡單的Gamma分布的密度函數(shù)：G a m m a ( t ∣ α ) = t α − 1 e − t Γ ( α ) Gamma(t|\alpha)=\frac{t^{\alpha-1}e^{-t}}{\Gamma(\alpha)}Gamma(t∣α)=Γ(α)tα−1e−t?

如果做一個變換t = β x t=\beta xt=βx，代入(1.3.1)式中，可得∫ 0 ∞ ( β x ) α − 1 e − β x d ( β x ) Γ ( α ) = 1 \int_0^{\infty} \frac{(\beta x)^{\alpha-1}e^{-\beta x}d(\beta x)}{\Gamma(\alpha)}=1∫0∞?Γ(α)(βx)α−1e−βxd(βx)?=1。

于是就得到Gamma分布的更一般形式：

G a m m a ( x ∣ α , β ) = β x α − 1 e − β x Γ ( α ) Gamma(x|\alpha,\beta)=\frac{\beta x^{\alpha-1}e^{-\beta x}}{\Gamma(\alpha)}Gamma(x∣α,β)=Γ(α)βxα−1e−βx?

其中， α 為Gamma分布的 shape parameter，主要決定了曲線的形狀;而 β為Gamma分布的 rate parameter 或者 inverse scale parameter，主要決定了曲線有多陡。Gamma分布的歸一化常數(shù)(使得Gamma值為1的點)恰為Γ \GammaΓ函數(shù)在點 α 處的值Γ ( α ) \Gamma(\alpha)Γ(α)。

Gamma分布的期望E = α β E=\frac{\alpha}{\beta}E=βα?、方差D = α β 2 D=\frac{\alpha}{\beta^2}D=β2α?。

下面看一下不同參數(shù)α , β \alpha,\betaα,β情況下的Gamma分布：

import numpy as np

from scipy.stats import gamma

from matplotlib import pyplot as plt

# Gamma分布(Gamma Distribution)

alpha_values = [1, 2, 3, 3, 3]

beta_values = [0.5, 0.5, 0.5, 1, 2]

color = ['b','r','g','y','m']

x = np.linspace(1E-6, 10, 1000)

fig, ax = plt.subplots(figsize=(12, 8))

for k, t, c in zip(alpha_values, beta_values, color):

dist = gamma(k, 0, t)

plt.plot(x, dist.pdf(x), c=c, label=r'$\alpha=%.1f,\beta=%.1f$' % (k, t))

plt.xlim(0, 10)

plt.ylim(0, 2)

plt.xlabel('$x$')

plt.ylabel(r'$p(x|\alpha,\beta)$')

plt.title('Gamma Distribution')

plt.legend(loc=0)

plt.show()

顯示結(jié)果如下：

同時，我們可以發(fā)現(xiàn)Gamma分布的概率密度和Poisson分布在數(shù)學(xué)上的形式具有高度的一致性。參數(shù)λ \lambdaλ的Poisson分布，概率為：P o s s i o n ( x = k ∣ λ ) = λ k k ! e − λ Possion(x=k|\lambda)=\frac{\lambda^k}{k!}e^{-\lambda}Possion(x=k∣λ)=k!λk?e−λ。

而在Gamma分布的密度函數(shù)中取α = k + 1 , β = 1 \alpha=k+1,\beta = 1α=k+1,β=1，可以得到：G a m m a ( x ∣ α = k + 1 , β = 1 ) = x k e − x k ! Gamma(x|\alpha=k+1,\beta=1)=\frac{x^ke^{-x}}{k!}Gamma(x∣α=k+1,β=1)=k!xke−x?。

可以看到這兩個分布在數(shù)學(xué)形式上是一致的，只是Poisson分布式離散的，Gamma分布式連續(xù)的，可以直觀認(rèn)為，Gamma分布式是Poisson分布在正實數(shù)集上連續(xù)化版本。

2. Beta 分布

2.1 首先要了解一下Beta 函數(shù)。

Beta 函數(shù)的定義：B ( α , β ) = ∫ 0 1 x α − 1 ( 1 − x ) β − 1 d x B(\alpha,\beta)=\int_{0}^{1}x^{\alpha-1}(1-x)^{\beta-1}\ dxB(α,β)=∫01?xα−1(1−x)β−1dx，其中α , β > 0 α,β>0α,β>0。

Beta函數(shù)性質(zhì)：

對稱性：B ( α , β ) = B ( β , α ) \Beta(\alpha,\beta)=\Beta(\beta,\alpha)B(α,β)=B(β,α)

與Γ函數(shù)的關(guān)系：B ( α , β ) = Γ ( α ) Γ ( β ) Γ ( α + β ) \Beta(\alpha,\beta)=\frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha+\beta)}B(α,β)=Γ(α+β)Γ(α)Γ(β)?，其中 α,β>0

2.2 Beta分布

對Beta函數(shù)做個變形，可以得到如下式子：∫ 0 1 x α − 1 ( 1 − x ) β − 1 B ( α , β ) d x = 1 \int_{0}^{1}\frac{x^{\alpha-1}(1-x)^{\beta-1}}{B(\alpha,\beta)}dx=1∫01?B(α,β)xα−1(1−x)β−1?dx=1

取等式左邊積分中的函數(shù)作為概率密度，就得到了Beta分布的密度函數(shù)：B e t a ( x ∣ α , β ) = x α − 1 ( 1 − x ) β − 1 B ( α , β ) Beta(x|\alpha,\beta)=\frac{x^{\alpha-1}(1-x)^{\beta-1}}{B(\alpha,\beta)}Beta(x∣α,β)=B(α,β)xα−1(1−x)β−1?。

可以發(fā)現(xiàn)Beta分布的歸一化常數(shù)(使得Beta值為1的點)恰為Beta函數(shù)在 (α,β) 處的值。

Beta 分布的期望E = α α + β E=\frac{\alpha}{\alpha+\beta}E=α+βα?、方差D = α β ( α + β ) 2 ( α + β + 1 ) D=\frac{\alpha\beta}{(\alpha+\beta)^2(\alpha+\beta+1)}D=(α+β)2(α+β+1)αβ?。

下面看一下不同參數(shù)α , β \alpha,\betaα,β情況下的Beta分布：

# Beta 分布(Beta distribution)

import numpy as np

from scipy.stats import beta

from matplotlib import pyplot as plt

alpha_values = [1/3,2/3,1,1,2,2,4,10,20]

beta_values = [1,2/3,3,1,1,6,4,30,20]

colors = ['tab:blue', 'tab:orange', 'tab:green', 'tab:red', 'tab:purple',

'tab:brown', 'tab:pink', 'tab:gray', 'tab:olive']

x = np.linspace(0, 1, 1002)[1:-1]

fig, ax = plt.subplots(figsize=(14,9))

for a, b, c in zip(alpha_values, beta_values, colors):

dist = beta(a, b)

plt.plot(x, dist.pdf(x), c=c,label=r'$\alpha=%.1f,\ \beta=%.1f$' % (a, b))

plt.xlim(0, 1)

plt.ylim(0, 6)

plt.xlabel('$x$')

plt.ylabel(r'$p(x|\alpha,\beta)$')

plt.title('Beta Distribution')

ax.annotate('Beta(1/3,1)', xy=(0.014, 5), xytext=(0.04, 5.2),

arrowprops=dict(facecolor='black', arrowstyle='-'))

ax.annotate('Beta(10,30)', xy=(0.276, 5), xytext=(0.3, 5.4),

arrowprops=dict(facecolor='black', arrowstyle='-'))

ax.annotate('Beta(20,20)', xy=(0.5, 5), xytext=(0.52, 5.4),

arrowprops=dict(facecolor='black', arrowstyle='-'))

ax.annotate('Beta(1,3)', xy=(0.06, 2.6), xytext=(0.07, 3.1),

arrowprops=dict(facecolor='black', arrowstyle='-'))

ax.annotate('Beta(2,6)', xy=(0.256, 2.41), xytext=(0.2, 3.1),

arrowprops=dict(facecolor='black', arrowstyle='-'))

ax.annotate('Beta(4,4)', xy=(0.53, 2.15), xytext=(0.45, 2.6),

arrowprops=dict(facecolor='black', arrowstyle='-'))

ax.annotate('Beta(1,1)', xy=(0.8, 1), xytext=(0.7, 2),

arrowprops=dict(facecolor='black', arrowstyle='-'))

ax.annotate('Beta(2,1)', xy=(0.9, 1.8), xytext=(0.75, 2.6),

arrowprops=dict(facecolor='black', arrowstyle='-'))

ax.annotate('Beta(2/3,2/3)', xy=(0.99, 2.4), xytext=(0.86, 2.8),

arrowprops=dict(facecolor='black', arrowstyle='-'))

#plt.legend(loc=0)

plt.show()

顯示結(jié)果如下：

注意到這里當(dāng)α = 1 , β = 1 \alpha=1,\beta=1α=1,β=1時，Beta分布退化成一個均勻分布。

此外，前面提到Beta分布的均值E = α α + β E=\frac{\alpha}{\alpha+\beta}E=α+βα?，也可以通過采樣的方法，在一個Beta分布中，采樣，計算均值。

這里我們?nèi)?alpha; = 2 , β = 1 \alpha=2,\beta=1α=2,β=1，理論上均值應(yīng)該是2 3 \frac{2}{3}32?，我們看一下采樣方法的結(jié)果：

import numpy as np

import numpy.random as nprnd

import scipy.stats as spstat

import scipy.special as ssp

import itertools as itt

import matplotlib.pyplot as plt

import pylab as pl

# Beta 分布的平均數(shù)

N = (np.arange(200) + 3) ** 2 * 20

betamean = np.zeros_like(N, dtype=np.float64)

for idx, i in enumerate(N):

betamean[idx] = np.mean(nprnd.beta(2, 1, i))

plt.plot(N, betamean, color='steelblue', lw=2)

plt.xscale('log')

plt.show()

print(spstat.beta(2, 1).mean(), 2.0 / (2 + 1))

結(jié)果顯示如下：

這里可以看到，隨著采樣點的增加，樣本點的均值也就更加的收斂，更加的接近2 3 \frac{2}{3}32?。 這樣，這個圖片的結(jié)果也符合大數(shù)定理，隨著采樣點的增加，只要樣本點無限大，那么最終的均值就會無限的接近它的期望值2 3 \frac{2}{3}32?.

Beta分布可以理解為它表示概率的概率分布—也就是在我們不知道一件事的概率是多少的時候，它能表示一個概率的所有可能值。

2.3 共軛分布

Beta分布其實是二項分布的共軛分布。

在貝葉斯概率理論中，如果后驗概率P(θ|x)和先驗概率P(θ)滿足同樣的分布律，那么，先驗分布和后驗分布被叫做共軛分布，同時，先驗分布叫做似然函數(shù)的共軛先驗分布。

貝葉斯參數(shù)估計的基本過程是：先驗分布+數(shù)據(jù)知識=后驗分布

因此可以得到：B e t a ( p ∣ k , n − k + 1 ) + B i n o m C o u n t ( m 1 , m 2 ) = B e t a ( p ∣ k + m 1 , n − k + 1 + m 2 ) Beta(p|k,n-k+1)+BinomCount(m1,m2)=Beta(p|k+m1,n-k+1+m2)Beta(p∣k,n−k+1)+BinomCount(m1,m2)=Beta(p∣k+m1,n−k+1+m2)。其中經(jīng)過了m次伯努利實驗，有 m1 個比 p 小， m2 個比 p 大。更一般的，對于非負(fù)實數(shù)α , β \alpha,\betaα,β，我們有如下關(guān)系：

B e t a ( p ∣ α , β ) + B i n o m C o u n t ( m 1 , m 2 ) = B e t a ( p ∣ α + m 1 , β + m 2 ) Beta(p|\alpha,\beta)+BinomCount(m1,m2)=Beta(p|\alpha+m1,\beta+m2)Beta(p∣α,β)+BinomCount(m1,m2)=Beta(p∣α+m1,β+m2)

以上式子實際上描述的就是Beta-Binomial共軛(Beta-二項分布共軛)。共軛意思是先驗和后驗都服從同一個分布形式。這種形式不變，我們能夠在先驗分布中賦予參數(shù)很明確的物理意義，這個物理意義可以延伸到后驗分布中進(jìn)行解釋，同時從先驗變換到后驗的過程中從數(shù)據(jù)中補(bǔ)充的知識也容易有物理解釋。(還有另一個好處是:每當(dāng)有新的觀測數(shù)據(jù)，就把上次的后驗概率作為先驗概率，乘以新數(shù)據(jù)的likelihood，然后就得到新的后驗概率，而不必用先驗概率乘以所有數(shù)據(jù)的likelihood得到后驗概率。)

2.4 Beta分布的應(yīng)用

那么我們簡單說個Beta-Binomial共軛的應(yīng)用。用一句話來說，beta分布可以看作一個概率的概率分布，當(dāng)你不知道一個東西的具體概率是多少時，它可以給出了所有概率出現(xiàn)的可能性大小。

舉一個簡單的例子，熟悉棒球運動的都知道有一個指標(biāo)就是棒球擊球率(batting average)，就是用一個運動員擊中的球數(shù)除以擊球的總數(shù)，我們一般認(rèn)為0.266是正常水平的擊球率，而如果擊球率高達(dá)0.3就被認(rèn)為是非常優(yōu)秀的。現(xiàn)在有一個棒球運動員，我們希望能夠預(yù)測他在這一賽季中的棒球擊球率是多少。傳統(tǒng)的頻率學(xué)派會直接計算棒球擊球率，用擊中的數(shù)除以擊球數(shù)，但是如果這個棒球運動員只打了一次，而且還命中了，那么他就擊球率就是100%了，這顯然是不合理的，因為根據(jù)棒球的歷史信息，我們知道這個擊球率應(yīng)該是0.215到0.36之間才對。對于這個問題，我們可以用一個二項分布表示(一系列成功或失敗)，一個最好的方法來表示這些經(jīng)驗(在統(tǒng)計中稱為先驗信息)就是用beta分布，這表示在我們沒有看到這個運動員打球之前，我們就有了一個大概的范圍。beta分布的定義域是 (0,1) 這就跟概率的范圍是一樣的。接下來我們將這些先驗信息轉(zhuǎn)換為beta分布的參數(shù)，我們知道一個擊球率應(yīng)該是平均0.27左右，而他的范圍是0.21到0.35，那么根據(jù)這個信息，我們可以取α = 81 , β = 219 \alpha=81,\beta=219α=81,β=219。(這樣取值可以從Beta的均值和分布考慮)。具體看下它的概率分布圖：

import numpy as np

from scipy.stats import beta

from matplotlib import pyplot as plt

# 擊球率為0.27左右時的Beta分布

x = np.linspace(0, 1, 1002)[1:-1]

fig, ax = plt.subplots(figsize=(10,6))

dist = beta(81, 219)

plt.plot(x, dist.pdf(x), c='b',label=r'$\alpha=%.1f,\ \beta=%.1f$' % (81, 219))

plt.xlim(0, .6)

plt.ylim(0, 16)

plt.xlabel('$x$')

plt.ylabel(r'$p(x|\alpha,\beta)$')

plt.title('Beta Distribution')

plt.legend(loc=0)

plt.show()

x 軸表示各個擊球率的取值， x 對應(yīng)的 y 值就是這個擊球率所對應(yīng)的概率。也就是說beta分布可以看作一個概率的概率分布。

那么有了先驗信息后，現(xiàn)在我們考慮一個運動員只打一次球，那么他現(xiàn)在的數(shù)據(jù)就是”1中;1擊”。這時候我們就可以更新我們的分布了，讓這個曲線做一些移動去適應(yīng)我們的新信息。因beta分布與二項分布是共軛先驗的(Conjugate_prior)。那么后驗是：B e t a ( α 0 + h i t s , β 0 + m i s s e s ) Beta(\alpha_0+hits,\beta_0+misses)Beta(α0?+hits,β0?+misses)。

其中α 0 \alpha_0α0?和β 0 \beta_0β0?是一開始的參數(shù)，在這里是81和219。所以在這一例子里，α \alphaα增加了1(擊中了一次)。β \betaβ沒有增加(沒有漏球)。這就是我們的新的beta分布B e t a ( 81 + 1 , 219 ) Beta(81+1,219)Beta(81+1,219)。我們看下?lián)糁幸淮魏笮碌膿羟蚵蔅eta分布：

import numpy as np

from scipy.stats import beta

from matplotlib import pyplot as plt

# 擊中一次后的擊球率Beta分布

x = np.linspace(0, 1, 1002)[1:-1]

fig, ax = plt.subplots(figsize=(10,6))

dist = beta(81, 219)

plt.plot(x, dist.pdf(x), c='b',label=r'$\alpha=%.1f,\ \beta=%.1f$' % (81, 219))

dist = beta(82, 219)

plt.plot(x, dist.pdf(x), c='r',label=r'$\alpha=%.1f,\ \beta=%.1f$' % (82, 219))

plt.xlim(0, .6)

plt.ylim(0, 16)

plt.xlabel('$x$')

plt.ylabel(r'$p(x|\alpha,\beta)$')

plt.title('Beta Distribution')

plt.legend(loc=0)

plt.show()

可以看到這個分布其實沒多大變化，這是因為只打了1次球并不能說明什么問題。但是如果我們得到了更多的數(shù)據(jù)，假設(shè)一共打了300次，其中擊中了100次，200次沒擊中，那么這一新分布就是：B e t a ( 81 + 100 , 219 + 200 ) Beta(81+100,219+200)Beta(81+100,219+200)。我們來看一下?lián)糁?00次，失誤200次之后新的擊球率Beta分布：

import numpy as np

from scipy.stats import beta

from matplotlib import pyplot as plt

# 擊中100，失誤200次之后的擊球率Beta分布

x = np.linspace(0, 1, 1002)[1:-1]

fig, ax = plt.subplots(figsize=(10,6))

dist = beta(81, 219)

plt.plot(x, dist.pdf(x), c='b',label=r'$\alpha=%.1f,\ \beta=%.1f$' % (81, 219))

dist = beta(181, 419)

plt.plot(x, dist.pdf(x), c='r',label=r'$\alpha=%.1f,\ \beta=%.1f$' % (181, 419))

plt.xlim(0, .6)

plt.ylim(0, 22)

plt.xlabel('$x$')

plt.ylabel(r'$p(x|\alpha,\beta)$')

plt.title('Beta Distribution')

plt.legend(loc=0)

plt.show()

注意到這個曲線變得更加尖，并且平移到了一個右邊的位置，表示比平均水平要高。

一個有趣的事情是，根據(jù)這個新的beta分布，我們可以得出他的數(shù)學(xué)期望為：α α + β = 0.303 \frac{\alpha}{\alpha+\beta}=0.303α+βα?=0.303，這一結(jié)果要比直接的估計要小100 100 + 200 = 0.333 \frac{100}{100+200}=0.333100+200100?=0.333。你可能已經(jīng)意識到，我們事實上就是在這個運動員在擊球之前可以理解為他已經(jīng)成功了81次，失敗了219次這樣一個先驗信息。

因此，對于一個我們不知道概率是什么，而又有一些合理的猜測時，beta分布能很好的作為一個表示概率的概率分布。

3. 狄利克雷分布(Dirichlet Distribution)

Dirichlet(狄利克雷)可以看做是將Beta分布推廣到多變量的情形。概率密度函數(shù)定義如下 :

D i r ( p ? ∣ α ? ) = 1 B ( α ? ) ∏ k = 1 K p k α k − 1 ( ∗ ) Dir(\vec p|\vec \alpha) = \frac{1}{B(\vec \alpha)} \prod_{k=1}^{K}p_{k}^{\alpha_{k}-1} \ \ \ \ \ (*)Dir(pundefined?∣αundefined)=B(αundefined)1?k=1∏K?pkαk?−1?(∗)

其中，α ? = ( α 1 , α 2 , … , α K ) \vec \alpha = (\alpha_{1},\alpha_{2},\ldots,\alpha_{K})αundefined=(α1?,α2?,…,αK?)為Dirichlet分布的參數(shù)。滿足：α 1 , α 2 , … , α K > 0 \alpha_{1},\alpha_{2},\dots,\alpha_{K} > 0α1?,α2?,…,αK?>0。

B ( α ? ) B(\vec \alpha )B(αundefined)表示 Dirichlet分布的歸一化常數(shù):

B ( α ? ) = ∫ ∏ k = 1 K p k α k − 1 d p ? B(\vec \alpha)=\int \prod_{k=1}^{K}p_{k}^{\alpha_{k}-1} \ d\vec pB(αundefined)=∫k=1∏K?pkαk?−1?dpundefined?

類似于Beta函數(shù)有以下等式成立：

B ( α ? ) = Γ ( ∑ k = 1 K α k ) ∏ k = 1 K Γ ( α k ) B(\vec\alpha) = \frac{\Gamma(\sum_{k=1}^{K}\alpha_{k})}{\prod_{k=1}^{K}\Gamma(\alpha_{k})}B(αundefined)=∏k=1K?Γ(αk?)Γ(∑k=1K?αk?)?

標(biāo)簽： gamma分布密度函數(shù) gamma函數(shù)可視化 gamma整數(shù)域表示 gamma實數(shù)域表示