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。
- Set b = e + k e ， e ≈ 2.71828 b = \frac{e + k}{e}，e \approx 2.71828 b=ee+k，e≈2.71828.
- 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.
-
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.
-
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.
- 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−u1u1).
- y = k e v y = ke^v y=kev.
- z = u 1 2 u 2 z = {u_1}^2u_2 z=u12u2.
- 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.
- 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
- b = ( e + 0.5 ) / e ≈ 1.184 b = (e + 0.5)/e \approx 1.184 b=(e+0.5)/e≈1.184.
- 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.
- 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.
- 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
- 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.
- 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.
- w + d − c z = − 1.47015 < 0 ， → w + d - cz = - 1.47015 < 0，\to w+d−cz=−1.47015<0，→ step 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.
- 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