伽马先验分布的草案

1先验分布服从伽马-逆伽马分布
1：取形状参数先验分布为伽马分布：
)exp()
(),(~)(1
bm m
a b
b a Ga m a a
-=
-Γπ，
其中a 为形状参数，b 为尺度参数。

2:尺度参数的先验分布为逆伽马分布：
⎪⎪⎭
⎫ ⎝⎛-
⎪⎪⎭
⎫
⎝⎛=+ηη
ηπv u v
v u IGa u u
exp 1
)(),(~)(1
Γ， 其中u 为形状参数，v 为尺度参数。

则两参数的联合验前分布为：
⎪⎪⎭
⎫ ⎝⎛-
⎪⎪⎭⎫ ⎝⎛⋅
-=
⋅=+-ηηηπv u v
bm m
a b
v u IGa b a Ga m u u
a a
exp 1)()exp()
(),(),(),(1
1
ΓΓ )(⋅Γ为伽马函数：
dx e x
a x
a -∞+-⎰
=
1
)(Γ 0>a
2 超参数确定方法
上面给出的先验分布中，除了两参数数都取无信息先验分布，其它先验分布都含有未知的超参数。

超参数可以根据经验专家给出，但是更多的时候要利用先验数据确定。

通常利用先验矩可以确定超参数。

这时，首先要获得参数θ(对于Weibull 而言，就是m 和η）的样本，然后才能去估计参数的样本矩。

然而，已知的先验数据与参数的样本并没有明显的对应关系。

在工程应用中，常用自助法（Bootstrap ）获得参数的样本。

自助法的核心是利用自助样本（或称为再生样本）来估计未知概率测度的某种统计量的统计特性。

设),,(21n t t t T ⋅⋅⋅=是得到的一组数控系统无故障工作时间样本。

通过其运用自助法便可获得m 和η的样本。

具体步骤如下：
1）对),,(21n t t t T ⋅⋅⋅=进行有放回抽样，可得到自助样本；
2）利用自助样本),,(*
*2*1*n t t t T ⋅⋅⋅=进行最大似然估计，得到m
ˆ，ηˆ； 3）重复上述两个步骤N 次，得到估计参数样本)}ˆ,ˆ(,),ˆ,ˆ(),ˆ,ˆ{(2211N N m m m
ηηη⋅⋅⋅； 4）利用得到的估计参数样本即可求未知参数m ,η的期望和方差。

通过得到的N 组参数样本，分别计算其期望与方差，可以得到验前分布的超参数。

当假设先验分布为二元正态时：
∑
==
N i i m m
N
1
ˆ1μ，∑
==N i i N
1
ˆ1ημη，∑=--=
N
i m i
m
m N
1
2
2
ˆ1
1）（μσ，
∑=--=
N
i i
N
1
2
2
ˆ1
1）（ηημησ，
当假设先验分布为伽马-逆伽马分布时，解下面两个方程组便可得。

⎪⎪⎩⎪
⎪⎨⎧--==∑∑==N
i m i N i i m N a
b m N a b
122
1
ˆ11ˆ1）（μ， ⎪⎪⎩⎪
⎪⎨⎧--=--=-∑∑==N
i i N i i N u u v N u v
122
21ˆ11)
2()1(ˆ11）（ημηη 3 Bayes 方法得到的后验分布
)
()
()|()|(t m t p t θπθθπ=
； （2.5） θ
θπθθ
d t p t M )()|()(⎰=
； （2.6）
由于)(t M 不依赖于θ，在计算θ的后验分布中仅起到一个正规化因子的作用。

假如把)(t m 省略，把Bayes 公式改写为如下等价的形式：
)()|()|(θπθθπt p t ∝
其中符号∝代表两边仅差一个常数因子，一个不依赖于θ的常数因子。

则联合后验分布为：
()()⎪⎪⎭
⎫ ⎝
⎛-
--⋅⎪⎪⎭
⎫ ⎝⎛⎪⎪⎭
⎫
⎝⎛⋅ΓΓ∝=
∑∏+∈-∈-+m
D C i i
D
i i b u r
m
u
b
t v am t m
m
u b v
a t M m m t L t m )(exp )
(1)()
,(),|()|,(1
1
1
ηηηη
ηπηηπβ （2.7）
()()dm
d t m t U L
U
L
m m ηηππηη
⎰⎰=
|,； （2.8）
在这里
()()
u b v
a u
b
ΓΓ是不依赖于βα,的常量因子，在以后的积分当中可以消掉，因此可以省略
它写成如下形式：
⎪⎪⎭
⎫ ⎝⎛---⋅⎪⎪⎭
⎫ ⎝⎛⎪⎪⎭
⎫ ⎝⎛=∑∏+∈-∈-+m D C i i m D
i i b u r
m
t v am t m
m
m H )(exp )
(1),(1
1
1
ηηηηη （2.9）
对上式进行积分可以得到：
()()dm
d m H J U L
U
L
m m ηηβαηη
⎰⎰=
,, （2.10）
4 参数m ,η和MTBF 值的点估计和区间估计 ()⎰
⎰=U
L
U
L
m m K
K dm d m H m
J
Z m E ηη
ηη,1)|( （2.11） ()J
J Z m E 1|=
; ()dm d m mH J U L
U
L
m m ηηηη
⎰⎰=
,1 （2.12）
m 的双侧区间估计⎪⎭
⎫
⎝⎛-
-m m ,，由以下公式求得：
()⎰⎰
-
=
U
L
L
m
m d m d m H J
ηηγ
ηη2
,1 （2.13）
()⎰⎰
-
-
=U
L
L
m
m d m d m H J
ηηγ
ηη2
1,1 （2.14）
同样的方法可以得到α的点估计和区间估计 ()⎰
⎰=
U
L
U
L
m m K
K
dm d m H J
Z E ηη
ηηη
η
,1)|( （2.15）
()J
J Z E 2|=
η ; ()dm d m H J U
L
U
L
m m ηηηηη
⎰⎰=
,2 （2.16）
η的双侧区间估计⎪⎭
⎫
⎝⎛--ηη,，由以下公式求得：
()⎰⎰-=
U
L
L
m m d m d m H J
η
η
γ
ηη2
,1 （2.17）
()⎰⎰-
-
=U
L
L
m m d m d m H J
η
η
γ
ηη2
1,1 （2.18）
MTBF 的点估计和区间估计
MTBF=t=⎪⎭
⎫
⎝
⎛
+
Γm 11η （2.19） ()⎰⎰=
U L
U
L
m m K
K
dm d m H t J
Z t E ηη
ηη,1)|( （2.20）
()3t |J E Z J
=
; ()dm d m tH J U L
U
L
m m ηηηη
⎰⎰=
,3 （2.21）
MTBF 的双侧区间估计，由以下公式求得：
()⎰⎰⎪
⎭⎫ ⎝
⎛
-Γ-
=
U
L
L
m m m t
dm d m H J
112
,1η
γ
ηη （2.22）
()⎰⎰⎪
⎭
⎫ ⎝⎛
-Γ-
-
=U
L
L
m m m t
dm d m H J
112
1,1η
γ
ηη （2.23）
5算例
通过新产品现场试验，我们得到一组现场样本T ：
完全失效数据
203
323
934
1583 1947 1961 3206
截尾数据 783 1353 1016 3097 3300
通过由自助法得到的m 和η样本，我们便可得到样本均值、样本方差：
m μ=1.44，ημ=2251，2
m
σ=0.025，2ησ=1.98×105
⎪⎪⎩⎪⎪⎨
⎧==025.044.12a
b a
b
， ⎪⎪⎩⎪
⎪⎨⎧⨯=--=-52
2
1098.1)
2()1(29621u u v u v
可得：95.82,6.57==b a ，134209,3.46==v u
m 的点估计（后验期望估计）
)(ˆm E m
E ==1.35 η的点估计（后验期望估计）
)(ˆηη
E E ==2888.3
取显著性水平α=0.1， 则m 的区间估计：
L m =1.16，U m =1.49
η的区间估计：
L η=2283.4，U η=3073.6
MTBF 值的点估计：
)ˆ11(ˆE E MTBF m
t +Γ=η
=2648。