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