gamma分布密度函數介紹（Gamma函數可視化什么意思）-環球速看

1. Gamma 分布

(資料圖)

1.1 首先要了解一下Gamma 函數。

Gamma 函數在實數域可以表示為：

Γ ( 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;

在整數域可以表示為：

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

此外，Gamma 函數具有如下性質：

Γ ( 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Γ函數是階乘在實數上的推廣。

1.2 Gamma函數可視化

"""

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

結果顯示如下：

當然這里只取了大于0的部分。

1.3 Gamma 分布

對Gamma函數做個變形，可以得到如下式子：∫ 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)

取等式左邊積分中的函數作為概率密度，就得到一個簡單的Gamma分布的密度函數：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分布的歸一化常數(使得Gamma值為1的點)恰為Γ \GammaΓ函數在點 α 處的值Γ ( α ) \Gamma(\alpha)Γ(α)。

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

下面看一下不同參數α , β \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()

顯示結果如下：

同時，我們可以發現Gamma分布的概率密度和Poisson分布在數學上的形式具有高度的一致性。參數λ \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分布的密度函數中取α = 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?。

可以看到這兩個分布在數學形式上是一致的，只是Poisson分布式離散的，Gamma分布式連續的，可以直觀認為，Gamma分布式是Poisson分布在正實數集上連續化版本。

2. Beta 分布

2.1 首先要了解一下Beta 函數。

Beta 函數的定義：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函數性質：

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

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

2.2 Beta分布

對Beta函數做個變形，可以得到如下式子：∫ 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

取等式左邊積分中的函數作為概率密度，就得到了Beta分布的密度函數：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?。

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

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

下面看一下不同參數α , β \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()

顯示結果如下：

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

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

這里我們取α = 2 , β = 1 \alpha=2,\beta=1α=2,β=1，理論上均值應該是2 3 \frac{2}{3}32?，我們看一下采樣方法的結果：

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 分布的平均數

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

結果顯示如下：

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

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

2.3 共軛分布

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

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

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

因此可以得到：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)。其中經過了m次伯努利實驗，有 m1 個比 p 小， m2 個比 p 大。更一般的，對于非負實數α , β \alpha,\betaα,β，我們有如下關系：

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-二項分布共軛)。共軛意思是先驗和后驗都服從同一個分布形式。這種形式不變，我們能夠在先驗分布中賦予參數很明確的物理意義，這個物理意義可以延伸到后驗分布中進行解釋，同時從先驗變換到后驗的過程中從數據中補充的知識也容易有物理解釋。(還有另一個好處是:每當有新的觀測數據，就把上次的后驗概率作為先驗概率，乘以新數據的likelihood，然后就得到新的后驗概率，而不必用先驗概率乘以所有數據的likelihood得到后驗概率。)

2.4 Beta分布的應用

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

舉一個簡單的例子，熟悉棒球運動的都知道有一個指標就是棒球擊球率(batting average)，就是用一個運動員擊中的球數除以擊球的總數，我們一般認為0.266是正常水平的擊球率，而如果擊球率高達0.3就被認為是非常優秀的。現在有一個棒球運動員，我們希望能夠預測他在這一賽季中的棒球擊球率是多少。傳統的頻率學派會直接計算棒球擊球率，用擊中的數除以擊球數，但是如果這個棒球運動員只打了一次，而且還命中了，那么他就擊球率就是100%了，這顯然是不合理的，因為根據棒球的歷史信息，我們知道這個擊球率應該是0.215到0.36之間才對。對于這個問題，我們可以用一個二項分布表示(一系列成功或失敗)，一個最好的方法來表示這些經驗(在統計中稱為先驗信息)就是用beta分布，這表示在我們沒有看到這個運動員打球之前，我們就有了一個大概的范圍。beta分布的定義域是 (0,1) 這就跟概率的范圍是一樣的。接下來我們將這些先驗信息轉換為beta分布的參數，我們知道一個擊球率應該是平均0.27左右，而他的范圍是0.21到0.35，那么根據這個信息，我們可以取α = 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 值就是這個擊球率所對應的概率。也就是說beta分布可以看作一個概率的概率分布。

那么有了先驗信息后，現在我們考慮一個運動員只打一次球，那么他現在的數據就是”1中;1擊”。這時候我們就可以更新我們的分布了，讓這個曲線做一些移動去適應我們的新信息。因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?是一開始的參數，在這里是81和219。所以在這一例子里，α \alphaα增加了1(擊中了一次)。β \betaβ沒有增加(沒有漏球)。這就是我們的新的beta分布B e t a ( 81 + 1 , 219 ) Beta(81+1,219)Beta(81+1,219)。我們看下擊中一次后新的擊球率Beta分布：

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次球并不能說明什么問題。但是如果我們得到了更多的數據，假設一共打了300次，其中擊中了100次，200次沒擊中，那么這一新分布就是：B e t a ( 81 + 100 , 219 + 200 ) Beta(81+100,219+200)Beta(81+100,219+200)。我們來看一下擊中100次，失誤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()

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

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

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

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

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

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分布的參數。滿足：α 1 , α 2 , … , α K > 0 \alpha_{1},\alpha_{2},\dots,\alpha_{K} > 0α1?,α2?,…,αK?>0。

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

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函數有以下等式成立：

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

標簽： gamma分布密度函數 gamma函數可視化 gamma整數域表示 gamma實數域表示