其他分享
首页 > 其他分享> > MCS:连续随机变量——Gamma分布

MCS:连续随机变量——Gamma分布

作者:互联网

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

    • 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).
    • 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.
    • 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.
    • 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.
  1. 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

在这里插入图片描述

标签:approx,frac,MCS,var,np,theta,随机变量,Gamma
来源: https://blog.csdn.net/u014281392/article/details/122847475