Gamma分布

    Gamma分布几乎跟Erlang分布一样，唯一的区别是参数 k k k，在Gamma分布中 k k k可以是大于0的任意值。

Gamma分布的概率密度函数:
f ( x ) = x k − 1 θ k e − θ x Γ ( k ) ， x > = 0 f(x) = \frac{x^{k-1}\theta^k e^{-\theta x}}{\Gamma(k)}，x >= 0 f(x)=Γ(k)xk−1θke−θx​，x>=0

Γ ( k ) \Gamma(k) Γ(k):Gamma函数:
Γ ( k ) = ∫ 0 ∞ t k − 1 e − t d t ， k > 0 \Gamma(k) = \int_0^{\infty} t^{k-1}e^{-t}dt， k > 0 Γ(k)=∫0∞​tk−1e−tdt，k>0
如果 k k k为正整数(Gamma = Erlang):
Γ ( k ) = ( k − 1 ) ! \Gamma(k) = (k-1)! Γ(k)=(k−1)!

Gamma变量的均值和方差:
μ = k θ \mu = \frac{k}{\theta} μ=θk​

σ 2 = k θ 2 \sigma^2 = \frac{k}{\theta^2} σ2=θ2k​

生成Gamma变量要比Erlang更复杂，要分两种不同的情况： k < 1 k < 1 k<1和 k > 1 k>1 k>1。

k < 1 k < 1 k<1(Accept-Reject Method)

x ′ x' x′: θ = 1 \theta = 1 θ=1的gamma变量， x x x：gamma变量， θ > 0 \theta > 0 θ>0。

1. Set b = e + k e ， e ≈ 2.71828 b = \frac{e + k}{e}，e \approx 2.71828 b=ee+k​，e≈2.71828.
2. Generate tow random uniform u ∼ U ( 0 , 1 ) ， u 1 , u 2 u \sim U(0,1)，u_1, u_2 u∼U(0,1)，u1​,u2​.
   • $ p = b \times u_1$
   • if p > 1 ， → p > 1，\to p>1，→ step 4.
   • if p < = 1 ， → p <=1，\to p<=1，→ step 3.
3. y = p 1 / k y = p^{1/k} y=p1/k.
   • if u 2 < = e − y u_2 <= e^{-y} u2​<=e−y，set x ′ = y x' = y x′=y， → \to → step 5.
   • if u 2 > e − y u_2 > e^{-y} u2​>e−y， → \to → step 2.
4. y = − l n ( b − p k ) y = -ln(\frac{b - p}{k}) y=−ln(kb−p​).
   • if u 2 < = y k − 1 u_2 <= y^{k-1} u2​<=yk−1，set x ′ = y x' = y x′=y, → \to → setp 5.
   • if u 2 > y k − 1 u_2 > y^{k-1} u2​>yk−1， → \to → 2.
5. Return x = x ′ θ x = \frac{x'}{\theta} x=θx′​

k > 1 k > 1 k>1

1.  
   • Set a = 1 2 k − 1 a = \frac{1}{\sqrt{2k - 1}} a=2k−1 ​1​.
   • b = k − l n 4 b = k - ln 4 b=k−ln4.
   • q = k + 1 / a q = k + 1/a q=k+1/a.
   • c = 4.5 c = 4.5 c=4.5.
   • d = 1 + l n ( 4.5 ) d = 1 + ln(4.5) d=1+ln(4.5).
2.  
   • Generate two random uniform u ∼ U ( 0 , 1 ) u \sim U(0,1) u∼U(0,1) variates, u 1 , u 2 u_1, u_2 u1​,u2​.
   • v = a × l n ( u 1 1 − u 1 ) v = a \times ln(\frac{u_1}{1 - u_1}) v=a×ln(1−u1​u1​​).
   • y = k e v y = ke^v y=kev.
   • z = u 1 2 u 2 z = {u_1}^2u_2 z=u1​2u2​.
   • w = b + q v − y w = b + qv - y w=b+qv−y.
3.  
   • if w + d − c z > = 0 w + d - cz >= 0 w+d−cz>=0，set x ′ = y x' = y x′=y， → \to → step 5.
   • if w + d − c z < 0 w + d - cz < 0 w+d−cz<0， → \to → step 4.
4.  
   • if w > = l n ( z ) w >= ln(z) w>=ln(z)，set x ′ = y x' = y x′=y， → \to → step 5.
   • if w < l n ( z ) w < ln(z) w<ln(z)， → \to → step 2.
5. Return x = x ′ θ x = \frac{x'}{\theta} x=θx′​.

例：随机变量 x x x服从伽马分布，其均值为：0.1，方差为：0.02。生成一个伽马变量。

μ = 0.1 = k θ ， σ 2 = 0.02 = k θ 2 \mu = 0.1 = \frac{k}{\theta}，\sigma^2 = 0.02 = \frac{k}{\theta^2} μ=0.1=θk​，σ2=0.02=θ2k​， θ = 5 ， k = 0.5 \theta = 5，k = 0.5 θ=5，k=0.5

1. b = ( e + 0.5 ) / e ≈ 1.184 b = (e + 0.5)/e \approx 1.184 b=(e+0.5)/e≈1.184.
2. u 1 = 0.71 ， u 2 = 0.21 ， p = 0.841 ， p < = 1 ， → u_1 = 0.71，u_2 = 0.21，p = 0.841，p <= 1，\to u1​=0.71，u2​=0.21，p=0.841，p<=1，→ step 3.
3. y = 0.707 ， e − y = 0.493 ， u 2 < = 0.493 ， x ′ = 0.707 ， → y = 0.707，e^{-y} = 0.493，u_2 <= 0.493，x' = 0.707，\to y=0.707，e−y=0.493，u2​<=0.493，x′=0.707，→ step 5.
4. x = 0.707 5 = 0.1414 x = \frac{0.707}{5} = 0.1414 x=50.707​=0.1414.

例：随机变量 x x x服从伽马分布，其均值为：10，方差为：66.6。生成一个伽马变量。

μ = 10 = k θ ， σ 2 = 66.6 = k θ 2 \mu = 10 = \frac{k}{\theta}，\sigma^2 = 66.6 = \frac{k}{\theta^2} μ=10=θk​，σ2=66.6=θ2k​， θ ≈ 0.15 ， k = 1.5 \theta \approx 0.15，k = 1.5 θ≈0.15，k=1.5

1. a ≈ 0.707 ， b ≈ 0.114 ， q = 2.914 ， c = 4.5 ， d ≈ 2.504 a \approx 0.707，b \approx 0.114， q = 2.914，c=4.5，d \approx 2.504 a≈0.707，b≈0.114，q=2.914，c=4.5，d≈2.504.
2. u 1 = 0.15 ， u 2 = 0.74 ， v ≈ − 1.226 ， y ≈ 0.440 ， z ≈ 0.0167 ， w ≈ − 3.899 u_1 = 0.15，u_2 = 0.74，v \approx -1.226，y \approx 0.440，z \approx 0.0167，w \approx -3.899 u1​=0.15，u2​=0.74，v≈−1.226，y≈0.440，z≈0.0167，w≈−3.899.
3. w + d − c z = − 1.47015 < 0 ， → w + d - cz = - 1.47015 < 0，\to w+d−cz=−1.47015<0，→ step 4.
4. − 3.899 > − 4.092 ， x ′ = 0.440 ， → -3.899 > -4.092，x'= 0.440，\to −3.899>−4.092，x′=0.440，→ step 5.
5. x = 0.44 0.15 ≈ 2.933 x = \frac{0.44}{0.15} \approx 2.933 x=0.150.44​≈2.933

Python模拟生成Gamma变量:

import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt

def generate_exponential_var(theta):
    u = np.random.uniform(0, 1)
    var = -1/theta*np.log(1 - u)
    return var

def generate_erlang_var(k=1, theta=1):
    if k == 1:
        var = generate_exponential_var(theta)
    else:
        var = 0
        for j in range(k):
            var += generate_exponential_var(theta)
    return var

def generate_gamma_var(k, theta):
    if k == 1:
        var = generate_erlang_var(k=1, theta=theta)
        return var
    elif k < 1:
        b = (np.e + k)/np.e
        while True:
            u1 = np.random.uniform(0, 1)
            u2 = np.random.uniform(0, 1)
            p = b*u1
            if p > 1:
                y = -np.log((b-p)/k)
                if u2 <= np.power(y, k-1):
                    x_ = y
                    break
                else:
                    continue
            else:
                y = np.power(p, 1/k)
                if u2 <= np.power(np.e, -y):
                    x_ = y
                    break
                else:
                    continue
    elif k > 1:
        a = 1/np.sqrt(2*k - 1)
        b = k - np.log(4)
        q = k + 1/a
        c = 4.5
        d = 1 + np.log(4.5)
        while True:
            u1 = np.random.uniform(0, 1)
            u2 = np.random.uniform(0, 1)
            v = a * np.log(u1/(1 - u1))
            y = k * np.power(np.e, v)
            z = np.power(u1, 2)*u2
            w = b + q*v -y
            if (w + d - c*z) >= 0:
                x_ = y
                break
            else:
                if w >= np.log(z):
                    x_ = y
                    break
                else:
                    continue
    var = x_/theta
    return var

μ = 1 ， σ 2 = 1 \mu = 1，\sigma^2 = 1 μ=1，σ2=1的Gamma变量： k = θ = 1 k = \theta = 1 k=θ=1 等价于 标准的指数分布:

gamma_2 : μ = 0.1 = k θ ， σ 2 = 0.02 = k θ 2 \mu = 0.1 = \frac{k}{\theta}，\sigma^2 = 0.02 = \frac{k}{\theta^2} μ=0.1=θk​，σ2=0.02=θ2k​， θ = 5 ， k = 0.5 \theta = 5，k = 0.5 θ=5，k=0.5

gamma_3 : μ = 10 = k θ ， σ 2 = 66.6 = k θ 2 \mu = 10 = \frac{k}{\theta}，\sigma^2 = 66.6 = \frac{k}{\theta^2} μ=10=θk​，σ2=66.6=θ2k​， θ ≈ 0.15 ， k = 1.5 \theta \approx 0.15，k = 1.5 θ≈0.15，k=1.5