基于 MCMC 的 invGamma 分布参数估计
摘要：文中主要介绍了马尔科夫链蒙特卡洛模拟在分布参数估计中的应
用。在理论部分，分析了 MCMC 等采样方法为什么可以从后验分布中
采样，而不能从似然函数中采样；在实证部分，对 inverse Gamma 分布
进行了参数估计，分别运用矩估计、最小二乘估计、极大似然估计和最
大后验估计方法进行估计，并运用 MCMC 方法从后验概率中采集样本，
用后验概率的均值再一次估计 invGamma 分布中的参数。
关键词： MCMC，逆 Gamma 分布，分布参数估计，
I
目录
0 引言 · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · 1
1 采样方法 · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · 1
1.1 似然函数中的采样问题 · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · 1
1.2 后验分布中的采样问题 · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · 2
1.3 Metropolis - Hasting 采样· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · 3
1.3.1 接受拒绝采样 · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · 3
1.3.2 Metropolis - Hasting 采样 · · · · · · · · · · · · · · · · · · · · · · · · · · · · · 4
2 基于 MCMC 的分布参数估计 · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · 5
2.1 invGamma 分布及构造样本· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · 5
2.2 矩估计法求分布参数 · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · 7
2.3 最小二乘法求分布参数 · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · 8
2.4 极大似然法求分布参数 · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · 9
2.5 最大后验概率法求分布参数· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · 10
2.6 MCMC 法从后验分布中采样 · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · 10
3 总结 · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · 12
参考文献 · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · 13
附录 · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · 17
II
0 引言
数学模型中一般会存在参数 (parameter)，并且参数估计问题一直是
各种模型的重要研究方向。考虑如下几种“方程式”的数学模型：用小写
x, y, · · · 表示变量，用大写的 X, Y, · · · 表示随机变量 (函数)，用 θ 表示模
型中得到参数
y = ax + b θ = (a, b) 方程 (组)
Y = ax + bX θ = (a, b) 随机变量方程 (组)
y
′
= ay + bx θ = (a, b) 常微分方程
a
∂
2
u
∂x
+ b
∂
2
u
∂y
+ c
∂
2
u
∂z
= 0 θ = (a, b, c) 偏微分方程
y
′
t
= ay
t−1
+ by
t−2
θ = (a, b) 时滞微分方程
Y
t
= aY
t−1
+ bY
t−2
+ cX
t
θ = (a, b, c) 离散随机微分方程
dX
t
= aX
t
dt + bX
t
dW
t
= (a, b) 随机微分方程
如果没有给出测量数据 (样本) 时，可以考虑给出一个具体的 θ，来对
模型进行模拟，以便观察模型。当给出具体的测量数据 D = {x
i
, y
i
}
n
i=1
时，
就要考虑利用数据 D 来估计模型中的参数 θ，一种易于想到的估计方法是
最小二乘法：min
θ
(y
i
− f (x
i
, θ))
2
。下面，我们将用一些估计方法来估计
分布 f 中的参数。值得一提的是，数据的分布函数 (或者密度函数) 本身就
是一个模型，分布参数的估计方法可以推广到上述模型当中，因此可以利
用分布参数估计问题来熟悉各种估计方法。在实验部分，我们将会利用矩
估计、最小二乘估计、极大似然估计、最大后验概率来估计 invGamma 分
布中的参数 θ = (α, β)，并用 MCMC 方法从最大后验概率中采样，用样本
均值继续估计参数。
1 采样方法
首先通过下面两个例子来说明为什么要采样。
1.1 似然函数中的采样问题
设样本数据为 D = {x
1
, · · · , x
n
}，并且分布为 f(θ )，x
i
iid
∼ f(θ)。如果采
用极大似然方法来估计参数 θ，要先构建似然函数 (样本的联合分布)，再
1
求 θ 使似然函数最大
L(θ, D) = f(x
1
, · · · , x
n
, θ) =
n
i=1
f(x
i
, θ)
其中：f(x
i
, θ) 是分布 f 在 x
i
处的概率，而 L(θ, D) 就是数据 D 的联合概
率分布 (参数 θ 的似然函数)。
我们的目标是找到 θ
∗
，使 L(θ) 最大。如果 θ 是一个随机变量，L(θ) 是
θ 的分布，那么求最大的过程就是求 θ 的众数，即 θ 的众数就是 θ
∗
。既然可
以使用 θ 的众数来进行估计，那么能不能考虑用 θ 的均值来做估计呢 (这
里仍然把 L(θ) 视为 θ 的分布)？如果均值不易求解，一个自然的想法是从
“分布”L(θ) 中采样，用样本均值代替总体均值进行估计。但要说明的一点
是，这里的 θ 是一个变量而不是随机变量，所以 L(θ) 并不是 θ 的分布，从
而并不一定能够从函数 L(θ) 采样，更谈不上均值。这也是为什么 L(θ, D)
称为 θ 的似然函数，称为 D 的联合概率分布的原因。
L(θ) 的函数示意图如图 (1) 所示
图 1: 似然函数示意图
这个地方，如果强行对 L(θ) 进行采样，应该仍然具有一定的可行性，
这个过程相当于用模拟退火算法来求 L(θ) 的最优解。
1.2 后验分布中的采样问题
在贝叶斯框架下，假设参数 θ 是随机变量，先验分布为 π(θ)，L(D|θ)
为条件分布 (本质是似然函数)，θ 的后验概率为
p(θ|D) =
π(θ)L(D|θ)
p(D)
=
π(θ)L(D|θ)
i
π(θ
i
)L(D|θ
i
)
2
其中：L(D|θ) ≜ L(D, θ) 就是数据 D 的联合概率分布，θ 的似然函数，π(θ)
为 θ 的先验分布。
和极大似然估计相似，我们的目标是求 θ 使得后验概率 p(θ|D) 最大。
由于
p(θ|D) =
π(θ)L(D|θ)
p(D)
∝ π(θ)L(D|θ)
从而
max
θ
p(θ|D) = max
θ
π(θ)L(D|θ)
求解上述最优化问题就是求后验分布 p(θ|D) 的众数。现在考虑用后验
分布的均值做估计，下面考虑如何计算复杂的后验分布的均值。由于 p(θ|D)
的复杂性，E =
∞
−∞
θp(θ|D)dθ 并不易于直接计算。如果能够求得后验分布
的样本 θ
1
, θ
2
, · · · 就好了，就可以用样本均值
¯
(θ) 来替代总体均值 Eθ，更
一般的，如果有了样本，不仅可以估计总体的分布，还可以估计总体的方
差、分布等等。因此，现在的问题就是如何从复杂分布 p(θ|D) 中采样。此
外，还要注意的是，在数据 D 给定下，每给一个参数 θ，L(D, θ) 是一个具
体的数值，π(θ) 也是一个具体的数值，但是后验概率 p(θ|D) 并能计算出来，
因为其分母 p(D) 是难于计算的。
1.3 Metropolis - Hasting 采样
1.3.1 接受拒绝采样
对于一个复杂的分布 f (θ)(这里要求分布 f 已知，每给定一个点 θ，都
会有一个具体的概率值)，用另一个易于采样的分布 g(θ) 的分布来生成新
样本 θ
∗
，然后判断是否接受 θ
∗
。这种想法对分布 g 有一些要求：
• g 要易于采样；
• g 的支撑集要包含 f 的支撑集。当产生不在 f 支撑集内的样本 θ
∗
时，
用概率 f(θ
∗
) = 0 来接受 θ
∗
；
• 如果 g 的形状与 f 的相似，则采样效率高，反之效率低下；
• f ≤ cg, ∀t。
3
下面考虑对于 g 生成的新样本 θ
∗
，应该以多大的概率接受呢？为方便，
这里假设 g 是均匀分布，c 为 f 的最大值，如图 (2) 所示。很明显，应该
以 f /cg 的概率接受 θ
∗
。这种做法的优点是，时间不随随机变量 θ 个数的
增加而增加，缺点是 g 不易于寻找。
图 2: 接受概率示意图
1.3.2 Metropolis - Hasting 采样
上面采样的关键问题是：如何产生新样本 θ，以及以多大的概率接受新
样本。如果一个随机过程 {θ
t
} 是 (宽/严) 平稳的，就可以将 {θ
t
} 的样本路
径视为平稳分布 π(θ
t
) 的一次采样，所以，如果一个随机过程 {θ
t
} 的平稳
分布是目标分布 f (θ)，那么样本路径即为 f(θ) 的样本。现在问题转化为如
何产生一个平稳分布为 f(θ) 的随机过程。
设置一个提议分布 q(·|·)(这里 q 并不是一个条件分布，而是将提议分
布 q 的参数写为条件的形式)，由提议分布产生新样本 θ
∗
t
，然后再以一定的
概率接受新样本，具体的 MH 算法如下 (1)
算法 1 Metropolis - Hasting for f(θ)
1: 初始化：初始状态 θ
0
，t := 1，t
max
。
2: for 对 t = 2, . . . , t
max
循环一下采样步骤 do
3: 产生新样本：θ
∗
t
∼ q(θ|θ
t
)；
4: 计算新样本的接受概率：α(θ
t
, θ
∗
t
) =
1,
f(θ
∗
t
)
f(θ
t
)
q(θ
t
|θ
∗
t
)
q(θ
∗
t
|θ
t
)
；
5: 判断是否接受新样本：如果 rand < α，则接受转移 θ
t
→ θ
∗
t
，即
θ
t+1
= θ
∗
t
，否则，不接受转移 θ
t+1
= θ
t
。
6: end for
上述算法对提议分布 q(·|·) 有一定的要求：
• q 要易于采样；
4
• q 的支撑集要包含 f 的支撑集；
• 给定一个具体的 θ
t
, θ
∗
t
，q 分布要能计算出具体的概率值 q(θ
t
|θ
∗
t
), q(θ
∗
t
|θ
t
)；
• 如果 q 的形状与 f 的相似，则采样效率高，反之效率低下；
• 在生成新样本 θ
∗
t
时，概率 q(·|θ
t
) 是依赖于 θ
t
的，而 θ
t
又是时变的，
所以要想采样效率高，要求 q 对其自身的分布参数不是很敏感。
另外，在计算过程中会出现 f(θ
t
), q(θ
t
) > 1 的情况，这是允许的并且注意
到接受概率 α 中包含
f(θ
∗
t
)
f(θ
t
)
，所以 MH 算法允许目标分布 f 中有可约的不
易计算的部分，比如后验分布 p(θ|D) =
π(θ)L(D|θ)
p(D)
中的 p(D) 项最终可以约
掉
p(θ
∗
t
|D)
p(θ
t
|D)
=
π(θ
∗
t
)L(D|θ
∗
t
)
π(θ
t
)L(D|θ
t
)
从而避免了复杂项 p(D)。
2 基于 MCMC 的分布参数估计
矩估计、最小二乘估计、极大似然估计、最大后验概率估计等参数估计
方法可以应用到各种模型当中，为了便于实验，下面用这些估计方法来估
计分布中的参数。以 invGamma 分布为例，因为没有实际的测量样本数据
D = {x
1
, x
2
, · · · , x
n
}，同时又为了便于观察估计量
ˆ
θ 是否接近真实 θ，设置
真实 θ = (α, β) = (3, 1)，然后从 invGamma(3,1) 中通过定义法产生样本数
据 D，再用 D 反过来估计分布参数。
2.1 invGamma 分布及构造样本
Gamma 分布的密度如下
f(x; α, β) =
1
Γ(α)β
α
x
α−1
e
−x/β
x > 0, α, β > 0
Gamma 分布的均值为 αβ，方差为 αβ
2
。设 X ∼ Gamma(α, β)，则 Y = 1/X
的分布定义为 invGamma 分布，由随机变量函数的密度计算公式可以计算
5
Y 的密度如下
f
Y
(y) = f
X
(1/y)
d
dy
y
−1
=
1
Γ(α)β
α
y
−α+1
exp(−1/βy)y
−2
=
(1/β)
α
Γ(α)
y
−α−1
exp(−(1/β)/y)
由此，invGamma(α, β) 的密度定义为
f(x; α, β) =
β
α
Γ(α)
x
−α−1
e
−
β
x
如果考虑从 invGamma(α, β) 中采样，考虑到 Gamma 分布和 invGamma
分布的关系:
X ∼ Gamma(α, β)
1
X
∼ inv Gamma (α, 1/β)
因此，如果 x ∼ Gamma(α, 1/β)，则 1/x ∼ invGamma(α, β)。所以只需要从
Gamma(α, 1/β) 中产生样本 x
i
，再去倒数 1/x
i
，这就是 invGamma(α, β)
中的样本。
通过上述方法生成 1000 个 invGamma(3, 1) 中的样本，样本频数如图
(3) 所示
图 3: 样本的频数直方图
6
2.2 矩估计法求分布参数
在有样本数据 D = {x
1
, x
2
, · · · , x
n
} 的情况下，使用矩估计来估计
invGamma(α, β) 中的参数。先来计算 invGamma 分布的均值和方差
EX = =
∞
0
xf(x; α, β)dx
=
β
α
Γ(α)
∞
0
x · x
−α−1
e
−
β
x
dx
=
β
α
Γ(α)
∞
0
x
−α
e
−
β
x
dx
=
β
α
Γ(α)
Γ(α − 1)
β
α−1
=
β
α − 1
(α > 1)
E
X
2
=
∞
0
x
2
f
(
x
;
α, β
)
d
x
=
β
α
Γ(α)
∞
0
x
−α+1
e
−
β
x
dx
=
β
2
(α − 1)(α − 2)
(α > 2)
当 α > 2 时，X 的方差为
DX = EX
2
− (EX)
2
=
β
2
(α − 1)
2
(α − 2)
联合 E, D 有
EX =
β
α − 1
DX =
β
2
(α − 1)
2
(α − 2)
用样本均值 ¯x 和方差 s
2
代替上述方程组的 E, D，解方程组后有
α
ME
=
(¯x)
2
s
2
+ 2
β
ME
=
(¯x)
2
s
2
+ 2
· ¯x
最终得估计结果为：(α, β)
ME
= (3.3378, 1.1419)。
7
2.3 最小二乘法求分布参数
最小二乘法常见于回归模型，下面将尝试使用最小儿乘法来估计分布
中的参数。将样本数据 D 最小到最大划分为 h 个区间 (h 为外来参数)，记
每个区间的中心为 x
i
，记落在第 i 个区间的样本频数为 n
i
，则第 i 个区间的
样本频率为 p
x
i
=
n
i
n
，再记区间大小为 ∆x
i
。如果仅考虑概率分布 f(x; α, β)
在 x
i
处的概率值 f(x
i
|θ)，则不易于求解
min
θ
h
i=1
[f(x
i
; θ) − p
x
i
]
2
这是由于样本频率 p
x
i
并不是用于近似 x
i
处的概率 f (x
i
)，而是用于近似
x
i
区间上的概率和 f(x
i
; θ)∆x
i
，所以最小二乘法的目标应该如下
θ
OLS
= arg min
θ
h
i=1
[f(x
i
; θ)∆x
i
− p
x
i
]
2
并且估计量 θ
OLS
与外来参数 h 有关，h 越大，估计越精确。OLS 估计得
结果为 (α, β)
OLS
= (2.8262, 0.9591)，估计结果如图 (4) 所示
图 4: 最小二乘估计结果
8
2.4 极大似然法求分布参数
用极大似然估计法估计 invGamma 分布中的参数。由于样本数据 D =
{x
1
, · · · , x
n
} 独立同分布，所以 x
1
, x
2
, · · · , x
n
的联合概率分布/似然函数为
L(D; θ) = P {x
1
, · · · , x
n
; θ}
=
n
i=1
f(x
i
; θ)
=
β
α
Γ(α)
n
n
i=1
x
−α−1
i
exp
n
i=1
−β
x
i
在 D 给定下，每给定一个参数 θ，都会有一个具体的数值 L(D; θ)。我
们的目标是求 θ 使得 L(D; θ) 最大
θ
MLE
= arg max
θ
L(D; θ)
= arg max
θ
ln L(D; θ)
上式中的 ln L(D; θ) 如下：
ln L = nα ln β − n ln Γ(α) +
n
i=1
(−α − 1) ln x
i
+
n
i=1
−
β
x
i
ln L(D; θ) 对 θ 求导，有
∂ ln L
∂β
=
nα
β
−
n
i=1
1
x
i
∂ ln L
∂α
= n ln β − n
Γ
′
(α)
Γ(α)
−
n
i=1
ln x
i
令上面的导数为 0，求解 α, β
nα
β
−
n
i=1
1
x
i
= 0
n ln β − n
Γ
′
(α)
Γ(α)
−
n
i=1
ln x
i
= 0
其中：
Γ
′
(α)
Γ(α)
= Φ(α) 为 Digamma 函数，其 MATLAB 命令为 psi，R 命令为
digamma。最终求解的极大似然估计结果为 (α, β)
MLE
= (2.8291, 0.9097)。
9
2.5 最大后验概率法求分布参数
假设参数 α, β 是相互独立的随机变量，并且先验分布 π(·) 为自由度为
2 和 3 的卡方分布
α ∼ χ
2
(2)
β ∼ χ
2
(3)
由贝叶斯公式，可以计算 α, β 的后验分布
p(θ|D) =
π(θ)L(D|θ)
p(D)
最大后验概率的目标是求 θ 使得后验分布概率最大，由于 p(θ|D) ∝ π(θ)L(D; θ)，
所以最优化模型变为
θ
MAP
= arg max
θ
p(θ|D)
=
arg max
θ
π(θ)L(D; θ)
最终求解的最优参数为 (α, β)
MAP
= (2.8222, 0.9075)。
上面的 θ
MAP
是后验分布 p(θ|D) 的众数，现在考虑用 p(θ|D) 的期望
来做 θ 的估计，下面用 MCMC 方法从后验分布 p(θ|D) 中采样。
2.6 MCMC 法从后验分布中采样
采用 Metropolis-Hasting 算法从后验分布 p(θ|D) 中采样，θ 共有 2 个
参数 α, β。设置提议分布为 χ
2
分布，并且逐个参数更新：如果 α 被更新
为 α
∗
，则在 α
∗
的基础上更新 β。Single-Component MH 算法如下 (2)
通过 MH 算法产生的样本路径如图 (5) 所示
10
算法 2 Single-Component Metropolis Hasting for p(θ|D)
1: 初始化：初始状态 θ
0
= [2, 0.5]，t := 1，马氏链长 T = 10000，θ
∗
t
= θ
t
= θ
0
，
马氏链记录 θ[1 : T ]。
2: for 对 t = 2, . . . , T 循环一下采样步骤 do
3: for i = 1 : 2 do // 对第 i 个参数而言
4: 产生新样本：θ
∗
t,i
∼ q(θ|θ
t,i
)；
5: 计算新样本的接受概率：α(θ
t
, θ
∗
t
) =
1,
f(θ
∗
t
)
f(θ
t
)
q(θ
t
|θ
∗
t
)
q(θ
∗
t
|θ
t
)
；
6: 判断是否接受新样本 θ
∗
t
：
7: if α(θ
t
, θ
∗
t
) > rand then
8: θ
t
= θ
∗
t
；
9: else
10: θ
∗
t
= θ
t
；
11: end if
12: end for
13: 更新马氏链：θ[t] = θ
t
。
14: end for
图 5: 样本路径图
从 t = 1000 之后的马氏链是收敛的，选取 1000 − 10000 作为后验分布
p(θ|D) 的样本。α, β 的后验分布如图 (6) 所示
11
图 6: mcmc 方法模拟的参数的后验分布
用后验分布的均值来估计分布参数，估计结果为 θ
P E
= (α, β)
P E
=
(2.8090, 0.9036)。
3 总结
不同方法估计得结果汇总如表 (1) 所示，上述方法都隶属于参数统计，
都是在总体分布已知的情况下进行参数估计 (总体分布已知，则样本的似然
函数就已知)。但是在许多复杂的模型中，似然函数和后验概率分布的最大
值是不容易求解的，或者说只能求解局部最值。因此，MCMC 方法相对而
言适应性更加广泛，并且 MCMC 采样后可以对目标分布有一个更好的掌
握，比如求目标分布的方差、分布等。但是 MCMC 方法有其自身的缺点，
就是相对而言计算量较大，需要多次模拟。
表 1: 估计结果汇总表
估计方法 α 的估计 β 的估计
真实值 3 1
矩估计 3.3378 1.1419
最小二乘估计 2.8262 0.9591
极大似然估计 2.8291 0.9097
最大后验概率估计 2.8222 0.9075
MCMC 后验均值估计 2.8090 0.9036
12
参考文献
[1] 杨爱军, 刘晓星, 林金官. 基于 MCMC 抽样的金融贝叶斯半参数 GARCH 模型研
究 [J]. 数理统计与管理.2015,34(03):452-462.
[2] 刘乐平, 高磊, 杨娜.MCMC 方法的发展与现代贝叶斯的复兴——纪念贝叶斯定理
发现 250 周年 [J]. 统计与信息论坛.2014,29(02):3-11.
[3] 鲁帆, 严登华. 基于广义极值分布和 Metropolis-Hastings 抽样算法的贝叶斯
MCMC 洪水频率分析方法 [J]. 水利学报.2013,44(08):942-949.
[4] 涂冬波, 蔡艳, 戴海琦, 丁树良.HO-DINA 模型的 MCMC 参数估计及模型性能研
究 [J]. 心理科学.2011,34(06):1476-1481.
[5] 卢斌, 华仁海. 基于 MCMC 方法的中国期货市场流动性研究 [J]. 管理科学学
报.2010,13(09):98-106+128.
[6] 赵琪.MCMC 方法研究 [D]. 山东大学.2007.
[7] 王权. 马尔可夫链蒙特卡洛”(MCMC) 方法在估计 IRT 模型参数中的应用 [J]. 考
试研究.2006(04):45-63.
13
附录
1 fu n c t i on main_MCMCmh
2 %% 人 为 设 计 的 分 布 参 数 为 alpha_real = 3 , beta = 1 。 下 面 绘 制 其 图 像
3 clc , c l e a r
4 alp ha_re al = 3 ; b eta_real =1;
5 xxx = 0 . 1 : 0 . 0 1 : 6 ;
6 yyy = invGamma_pdf( xxx , alpha_real , be ta_real ) ;
7 f i g u r e (1 )
8 p l o t ( xxx , yyy ) % 0 . 5 , 0 . 6 , 0 . 7 处 取 得 最 大 值
9 t i t l e ’ 实 际 密 度 ’
10
11 %% 人 为 构 建 样 本 sample
12 % 定 义 法 生 成invGamma分 布 的 随 机 数
13 T = 1000; % 采 样 数
14 x = gamrnd( alpha_real , 1 . / beta_real , 1 , T) ;
15 sample_invGam_define =1./x ;
16 f i g u r e (2 )
17 h i s t ( sample_invGam_define , 200)
18 t i t l e ( ’ 定 义 法 采 集 的invGamma分 布 的 样 本 ’ )
19 c l e a r x
20
21 %% 分 布 参 数 alpha 、 beta 的 矩 估 计
22 % sample_define 就 是 样 本
23 alpha_ME = mean( sample_invGam_define ) ^2/ var ( sample_invGam_define ) +
2 ;
24 beta_ME = (mean( sample_invGam_define ) ^2/ var ( sample_invGam_define )
+ 1) * mean( sample_invGam_define ) ;
25
26 %% 分 布 参 数 alpha 、 beta 的 最 小 二 乘 估 计 o l s
27 % 计 算 x i 以 及 x i 处 的 频 数
28 h = 200;
29 [ count_xi , x i ] = h i s t ( sample_invGam_define , h) ;
30 Pxi = count_xi/T; % xi 点 处 的 频 率 值
31
32 delta _xi = (max( sample_invGam_define )−min ( sample_invGam_define ) ) /h ;
33
34 % o l s 最 优 化
35 fun_ols = @( t heta ) sum( ( invGamma_pdf( xi , the t a (1) , theta (2 ) ) *
del ta_x i − Pxi ) . ^ 2 ) ; % 最 小 二 乘 目 标 ： sum( Fxi−Pxi ) ^2
36 the ta0 = [ 2 , 0 . 5 ] ; % 初 始 值
37 theta_OLS = fm ins earch ( fun_ols , t heta0 ) ;
38
14
39 alpha_OLS = theta_OLS (1) ;
40 beta_OLS = theta_OLS ( 2 ) ;
41
42 f i g u r e (5 )
43 p l o t ( xi , Pxi , ’ ro ’ )
44 hold on
45 p l o t ( xi , invGamma_pdf( xi , alpha_real , bet a_rea l ) * delta_xi , ’b− ’ )
46 p l o t ( xi , invGamma_pdf( xi , alpha_OLS , beta_OLS ) * delta_xi , ’ g<’ )
47 hold o f f
48 leg end ( ’ 样 本 频 率 ’ , ’ 实 际 概 率 ’ , ’ 估 计 概 率 ’ )
49 t i t l e ( ’ 最 小 二 乘 求 参 数 ’ )
50
51 c l e a r x i count_xi
52 %% 分 布 参 数 alpha 、 beta 的 极 大 似 然 估 计MLE
53 fun_mle = @( t heta ) MLE_eq_invGamma( theta , sample_invGam_define ) ;
54 the ta0 = [ 2 , 0 . 5 ] ;
55 theta_MLE = f s o l v e ( fun_mle , t het a0 ) ;
56
57 alpha_MLE = theta_MLE ( 1 ) ;
58 beta_MLE = theta_MLE ( 2 ) ;
59
60 %% 分 布 参 数 alpha 、 beta 的 最 大 后 验 估 计MAP
61 % the ta 的 先 验 pri or_pdf
62 hyperparameter = [ 2 , 3 ] ;
63 p r i o r = @( t h eta ) prior_pdf ( theta , hyperparameter ) ; %返 回 thet a 处 的 先
验 概 率 值
64
65 % 似 然 函 数 l i k e l i h o o d
66 li k e l i ho o d _ f u n = @( t h eta ) l i k e l i h o o d ( sample_invGam_define , theta ,
@invGamma_pdf) ;
67
68 % 后 验 概 率 p o s t e r i o r 即 为 目 标 分 布
69 p o s t erior_pd f = @( t heta ) p r i o r ( t heta ) * li k e l i ho o d _ f u n ( thet a ) ;
70
71 % MAP（极 大 换 为 极 小 ）
72 the ta0 = [ 2 , 0 . 5 ] ; % 初 始 值
73 obj_pdf = @( the ta ) − p r i o r ( t heta ) * l i k e l i ho o d _ f u n ( thet a ) ;
74 theta_MAP = fminsear ch ( obj_pdf , t heta 0 ) ;
75
76 alpha_MAP = theta_MAP( 1) ;
77 beta_MAP = theta_MAP(2 ) ;
78
79 %% mcmc参 数 估 计 后 验 的 均 值 或 众 数
15
80 T = 1000 0; % 采 样 数
81 n_theta = lengt h ( thet a0 ) ;
82 theta_mcmc = z er o s ( n_theta , T) ;
83 t = 1 ;
84 theta_mcmc ( : , t ) = theta0 ’ ;
85 xt = theta0 ’ ;
86 xt_new = xt ;
87
88 % 提 议 分 布 （ 匿 名 函 数 ） 卡 方 分 布
89 prop_pdf = @(xt_new , xt ) chi2pd f ( xt_new , xt ) ;
90
91 f o r t = 2: T
92 f o r i = 1 : n_theta
93 % 从 提 议 分 布 中 产 生 一 个 新 样 本 xt_new
94 xt_new ( i ) = prop_pdf_rnd ( xt ( i ) ) ;
95
96 % 计 算 新 样 本 xt_new的 接 受 概 率 accept_prob
97 accept_prob = min ( [ 1 , p o s t e rior_pdf (xt_new) * prop_pdf
( xt ( i ) , xt_new( i ) ) . . .
98 / ( posterio r _ p d f ( xt )
* prop_pdf ( xt_new( i ) , xt ( i ) ) ) ] ) ;
99
100 % 以 概 率 alpha_new 接 受 新 样 本 xt_new
101 i f accept_prob > rand
102 xt = xt_new ;
103 e l s e
104 xt_new = xt ;
105 end
106 end
107 theta_mcmc ( : , t ) = xt ;
108 end
109
110 c l e a r i t xt xt_new
111
112 % 判 断 样 本 是 否 收 敛
113 f i g u r e (3)
114 p l o t ( 1 :T, theta_mcmc )
115 leg end ( ’ alpha ’ , ’ beta ’ )
116 t i t l e ( ’ 样 本 路 径 tac e pl o t ’ )
117
118 % 从 1000 个 路 径 之 后 设 置 为 样 本
119 theta_sample = theta_mcmc ( : , 1001 : end ) ;
120
16
121 % 自 相 关 性 （ 相 关 性 很 强 ）
122 f i g u r e (4 )
123 subpl ot 211
124 a u tocorr ( theta_sample ( 1 , : ) )
125 subpl ot 212
126 a u tocorr ( theta_sample ( 2 , : ) )
127
128 % the ta 的 后 验 分 布 （ 频 率 直 方 图 ）
129 [ count_alpha , a l p h a i ] = h i s t ( theta_sample ( 1 , : ) , 10) ;
130 [ count_beta , b e t a i ] = h i s t ( theta_sample ( 2 , : ) , 10) ;
131 P_alpha = count_alpha /9000;
132 P_beta = count_beta /9 000;
133 f i g u r e (5)
134 subpl ot 121
135 bar ( alph ai , P_alpha )
136 t i t l e ’ alpha 的 后 验 分 布 ’
137 subpl ot 122
138 bar ( b etai , P_beta)
139 t i t l e ’ beta 的 后 验 分 布 ’
140 c l e a r count_alpha count_beta a l p h ai b e t a i
141
142 % the ta 后 验 分 布 的 均 值
143 alpha_mcmc = mean( theta_sample ( 1 , : ) ) ;
144 beta_mcmc = mean( theta_sample ( 2 , : ) ) ;
145
146 save ( ’ result_mcmc . mat ’ )
147
148
17
