229
龚光鲁，钱敏平著 应用随机过程教程及其在 算法与智能计算中的应用
清华大学出版社,2003
第 9 章 以图象信息为背景的随机场,迭代 Markov 系统
随机场的一个直观背景,是以图像为基本事件的随机元素,另一方面,迭代 Markov 系统是取值为压缩映射 随机映射 的迭代,其目的是得到其 不变元素,以用来描述一个图像,
1 有限格点上的 Markov 随机场与图象
1,1 有限格点上的 Markov 随机场
定义９,１ （ 相邻系统与子团 ）
设 G 为 平面有限格点 集合,}1,1:),{( 21 NjNijiG ≤≤≤≤=,对于 G 中的任意两点 yx ≠,可以定义它们相邻,或者不相邻,这样定义了 S 中的一种相邻概念以后,就得到了 G 上的一个 相邻系统,注意,两点是否相邻完全可以是抽象的,并不需要在直观上比邻
( 有时为了看起来有对称性,也可以认为 1N 与 1 恒同,2N 也与 1 恒同,即认为 G 为环面格点,)
在本书中,我们把 yx,相邻记为 yx?~,记 与 x 的相邻 的 点的 全体 为 )(x?,G 的子集
C 称为 G 的 子团,如果 C 中的任意两个元素都相邻,含 x 的子团 记为
U )()( xxC?=,
一个相邻系统由 G 中点对的相邻概念所决定,因此 也由全体子团所决定,所以人们常用? 来代表一个相邻系统,或用全体子团 }{C 来代表一个相邻系统,
例 ９．２ (1) 若 f=? )(x,则只有单点集是 G 的 子团,
(2) 若 }{\)( xGx =?,则 G 的所有子集都是 子团,,
(3) 一阶相邻 ( 最紧相邻,简称 紧邻 ),? 的相邻点 o可以 如下图
o
oo
o
此时 G 的全部 子团 形式 为,
f,?,,
|
(4) 二阶相邻,? 的相邻点 o可以 如下图
230
ooo
oo
ooo
此时 G 的全部 子团 形式 为,
f,?,,
|,
\ 及其转动,
\| 及其转动,
×
||,
高 阶相邻等 可以 如 下 图 所 示
6
54345
42124
631136
42124
54345
6
其中数字代表分阶 为 相邻的阶,S 的一个相邻系统也就是以 G 的点为顶点的一个双向图,
定义９,３ （ 随机场 ）
以 平面矩形格点集 G 中的格点 为参数,取值于 }1,,0{?LL ( 也可以是 实直线 R,或一个离散的集合 ) 的随机变量族,称为 取值于 }1,,0{?LL 的 G － 随机场,简称为 随机场,
按照 统计物理中的术语,我们称 }1,,0{?LL 为 自旋空间,
于是,随机场的一个样本,就是在格点集 G 的每个格点上,坐着,}1,,0{?LL 中某 一个 数,随机场的样本的一个直观 背景是图像．黑白图象是连续变化的一个灰度图,而 在实际图象处理中的 黑白灰度都是经过离散化采样的,设黑白灰度可分成 L 个 (一般 L=256)不同的等级,}1,,0{?LL 中的一个数值可以解释为黑白的灰度,最粗放的是黑白轮廓图，在每个采样点上只有,黑,,“白,两个象素之一 (即 1=L 的情形 ),一般 G 为 256× 256 个采样点 (即 25621 == NN ),那么,一幅黑白灰度图就可以看成在 G 的每一个格点上的一个灰度 ( 象素,Pixel),即给每一个格点赋予 1,,0?LL 中的一个值,所以,G － 随机场的样本,
即做一次试验可能得到的结果,可以解释为 G 上的一个黑白灰度图,
在这个意义上,随机场就也可以解释为 随机图,或者说,它是在图像空间取值的一个随机元素 ( 简称为 随机元 ),这是随机场的第 2 种解释,从相对灰度 ( 每个格点上的灰度除以总灰度 ( 灰度的和 )) 的角度也可以说,随机场的分布是格点集 G 上的一个随机的概率分布,
定义９,４ （组态） 从数学的角度抽象,随机场的一个样本就是 GLS }1,,0{?=? L 的一个元素,记为 Lba,,称为一个 组态 ( configuration ),即 ):( Gxx ∈= aa,
231
}1,,0{?∈ Lx La 称为组态 a 的 x 坐标,再记 ))(:( GJxxJ?∈= aa,称为组态 a 的 J
坐标,
所以,随机场也可以看成 是一个 随机组态,或者看成 为在 组态空间 GLS }1,,0{?= L
上取值的随机元,这是随机场的第 3 种解释,从分布的角度 也可以说,随机场 的 分布 是 组 态空间上的一个概率分布,
由此可见,一个灰度图就可以 近似地 视为 一个组态,
定义９,５ （ Markov 随机场 ）
随机场 }:{ Gxx ∈x 称为 Markov 场,如果对于任意格点 Gyx ∈,及任意组态
GLS },,1{ L=∈a,恒有
),|( xyP yyxx ≠== 一切axax )~,|( xyP yyxx
=== 一切axax,(9,1)
Markov 场 的含义是：在 任意格点 x 的余格点位置 上取 已知 值的条件下，随机场在格点 x 处取值的概率规律只与格点 x 的? 相邻点有关,
由定义 可以 直接 证明 下述定理
定理９．６ 随机场 }:{ Gxx ∈x 是 Markov 场 的 等价 条件为,对 任意格点集
GA? 及任意组态 GLS }1,,0{?=∈ La,恒有
),|,( AyAxP yyxx?=∈= axax ))(,|,( AyAxP yyxx?∈=∈== axax,
其中 }{)( 中点的邻点全体AA?=?,
1,２ 相邻系统 的 Gibbs 分布与 Gibbs 随机场 ( 邻位势 Gibbs 场 )
定义９．７ （ 格点集 G 的 子集与组态的相互作用 ）
对于 给定 的 一个子 格点集 )( GA?,组态上的一个 与之相关的 函数 )(aAU
( GL }1,,0{?∈ La ) 称为 A 对于组态 a 的相互作用,如果 它 满足,对于任意 两个组态
GL }1,,0{,?∈ Lba,只要 )()( bxax
xx = ( Ax ∈? ),就有 )()( ba AA UU = ( )(aAU
代表组态 a 在 A处的,坐标,对 a 的作用,所以 要求对在 A 处,坐标,相同 的 任意 组态的相互作用 一 样 ),
定义９．８ （ 位势函数 ） 在补充定义 0=fU 后,相互作用
}:{ GAU A?
的 全体 称为 位势,对于给定的相邻系统?,如果对于任意 A，只要它不是 子团，就有
0=AU,那么 }:{ 子团的为GCUC 就 称为 邻位势,如果 相邻正好就是直观上的 紧邻,则这样的 的 邻位势 简称 紧邻位势,
232
定义９,９ （ 能量函数 ） 对于 邻位势,组态的函数
∑
=
C
CU UH
子团
(9,2)
称为 邻位势的能量函数,
定义９,１０ （ 邻位势 的 Gibbs 分布与 Gibbs 随机场 ）
在组态空间 GLS }1,,0{?= L 上定义的概率分布 p ∈= aap (),( S ),称为 邻位势 ( 在能量函数下 ) 的 Gibbs 分布 ( 这里是离散分布 ),其中
)(1)( aap UHe
Z
=,
UH 为 邻位势 的能量函数,而 Z 为一个归一化常数,称为 Gibbs 分布 的 配分常数,
∑
∈
=
S
HUeZ
a
a )(,(9,3)
以 Gibbs 分布 为分布 的随机场 称为 Gibbs 场,即
)(),( apax =∈?= GxP xx,
就是说,Gibbs 分布 p是随机场在所有的格点处取值的联合分布,
在 第 6,7 章中的 Ising 模型中 的 Gibbs 分布,就是这里的 Gibbs 场的 概率分布的 特殊情形,即 2=L,自旋空间为 }1,1{? 的情形,其 对应的位势是紧邻的 ( 而把 外场看成格点对自己的作用的和 ),Ising 模型中 Gibbs 分布中有一 个倒温度参数,这是为了调节用的,这里无妨 省 去,
再则,由于在不同的邻域中 的科学家有 他们习惯用的记号,所以,在本节中的记号与第
6,7 章中用的有些不同,请读者见谅由 于 记号改动所带来的不便,我们列出其对比,
pba
bpahx
,
,,7,6
本节章第倒温度配称分布组态
t
tX
x
随机场
,
在第 6,7 章中,由于 Ising 模型较为简单,我们只提出了 Gibbs 分布 的概念,而并未提出随机场 的 概念,事实上,Gibbs 分布是 Gibbs 随机场的 不变 分布 （ 有的文献中 Gibbs 分布与
Gibbs 随机场的术语是混用的）,因此 在第 6,7 章中,我们如很多书上那样用 hx,表示组态,
而在本节中,我们强调的是随机场,它用 hx,表示更为合适,所以我们就与很多书上不一样,
改用 ba,表示组态,
定理 9.1 １ (Hammersley - Clifford 基本 定理 )
对于给定的相邻系统?, Markov 随机场与 邻位势 Gibbs 场 是 一样的,】
这个定理的意义在于,用 Markov 随机场来描述图象是一个直观认识,而用 邻位势 Gibbs 场 则 可 归结为其能量函数,后者有时可以用简单的参数来表达 ( 例如相互作用具有平移不变性的时候 ),那时就可用几个参数来描述图象,从而可以达到信息压缩的目的,于
233
是,此定理将直观认识转化为参数表示,以便于数学处理,
１,３ 图象处理的随机过程方法的思路 原则 概述
正如前面所述，一幅黑白图象 是 某个组态空间中的一个组态,在图像传输过程中,需要对于灰度图的信息进行处理与压缩,对于接收到的信息还需要恢复（ Retrieving ）为原图像．常见的方法有,
第 1 种方法 ( 用能量函数的 Gibbs 分布对图像建模 )
按黑白图像的相对灰度 ( 即每个格点的灰度除以图像的总灰度和 )大小,把图像看成格点 集 G 上的组态空间 GLS },,1,0{ L= 中的一个组态，或者看成 G 上的一个概率分布,由于自然景物等的图象往往具有一些局部相联系 ( 空间不变性 ) 的结构,这些局部结构可以抽象地归结为组态间按某种相邻系统的能量函数,
从而可以把一个图像用一个能量函数取到最小 值 的组态，也就是用 Gibbs 分布 达到 最大 值 处的组态，从而可以通过能量函数或 Gibbs 分布来对图像建模,这个 方法 的要害处是：如何 由 图像给出其能量函数,能 量函数表示图像能起一定的信息压缩作用．而反过来由能量函数恢复图像的问题,则 可以用模拟退火或其它
Markov Monte Carlo 方法运作,
第 2 种方法 ( 用低阶 Markov 场的样本对图像建模 )
将黑白图像看成随机场的一个样本,由定理９,１１知道任意随机场都可以看成 Markov 场,而图像的局部结构又说明,用较低阶相邻的 Markov 场的样本来描述图象更为合理,
第 3 种方法 (用迭代随机映射的不动点对图像建模 )
对于一个黑白轮廓图像，构造 与之相系的 有限个压缩映射，使得由此 压缩映射组 给出 的 一个迭 代映射的不动点近似地是此轮廓图．而对于一个黑白灰度图，构造与之相系的一个随机地取有限个压缩映射的随机映射，使此图像近似地是此 随机 映射的不变分布．（ 这 在第３节中将较多地介绍）,这种方法实际上是处理图像的分形方法的另一角度的发展,
一般地，彩色图可以看成红，黄，兰三张灰度图的叠置,
收到了不清晰的图像,也就是 受 到了干扰的图像,就需要对图像进行滤波,称为 图像的清洗,其中常见的一种清洗方法为 Bayes 方法,即把未清洗的图像看成先验分布,而寻找转移机制,使清洗后的恢复图象是此先验分布的 Bayes 分布 (后验分布 ),图像受到污染的典型情形,例如,被 测量的图象象素 xa 受到检测与录制系统的干扰,使 实际 测到的象素为
∑
∈
+=
)(
)(
xy
xxyx whg ab,
其中 g 为一个已知函数,yh 为已知的通道系数,而噪声 xw 的统计规律为已知,这时要恢复图象象素 xa
就需要进行滤波,
此外,图象边 缘信息也可能受到干扰而产生变形,也需要用边缘检测来恢复图象的信息,因此有时应在图象的先验信息中加进,边 -随机场,的先验分布,
有时 在处理图像前可能 还需要对图像进行 有部分相重叠的 分割 (segmentation),
１,４ Gibbs 分布的 样本的 Gibbs 采样法 (Gibbs Sampler)
要对 Gibbs 场作 一些统计 估计 （例如参数估计）,就 先 要得到 Gibbs 分布的配分常数,但是由于组态空间 GLS },,1{ L= 的元素数通常非常 庞大,例如 可以 有 61010,这就使 直接 计算
234
配分常数 Z 时 的求和项太多，从而 常会出现大量被求和的项的值小于计算机误差的精度，
以致使计算 在实际上 失去 可能,用 Markov 链 Monte Carlo 方法 (M C MC) 可以解决这个困难，
即先 近似地得到 Gibbs 分布,再 进行参数 估计,所以 MCMC 是 一个重要的计算手段,
M C MC 方法就是构造以 Gibbs 分布为可逆分布的 Markov 链的转移矩阵,因为这里的随机场的 Gibbs 分布远比 Ising 模型的 Gibbs 分布 复杂,我们只构造离散时间的 Markov 转移矩阵,而且还希望如 Glauber 动力学那样,每次只在一个格点上改变 （ 称为 按 Gibbs 方式 转移 ）,
最后还要保证由这样得到的转移矩阵的各行都收敛于 Gibbs 分布,
为了使构造更为透明,我们先对任给的格点子集 J 定义 如下的,部分转移,
=?
)(0
),()(1),( )( \
其它外坐标相同在 JeZp JSJUH
JJ
baaba ab,(9,4)
其中 JSJ \ab 表示,把组态 a 中在 J 处的坐标换成组态 b 中的相应坐标,而 )(aJZ 是一个归一化常数 (,局部 的,配分常数 ),
∑?=
J
JSJUH
J eZ
b
aba )( \)(,(9,5)
这种,部分转移,表示 组态间只在 位置 J 外的 坐标 相同时才可能 转移,
在此基础上 我们 定义转移矩阵 如下,记 G 中的全部格点的 一个给定的排序为
)|(|),,,{ 21||1 NNGxx G =L,任给组态 S∈ba,,定义
P J SJp ∈= baba,)),((,=P ∏
=
||
1
G
k
P }(
kx
,(9,6)
则 P是一个转移矩阵,注意由于 P }(
kx
中的转移只在一个格点 ix 上进行,由 P构造的 Markov
链 的一次转移,可以视为 || G 次连续转移的结果,其中 每次转移只在一个格点上进行,而且要做遍所有 的格点,
下面的 遍历 收敛 定理是 Ma rkov 链 的 Monte Carlo 方法的 理论 基础,而 Markov 链 Monte
Carlo 方法主要用于获得 Gibbs 分布 的样本,以便进一步得到 Gibbs 场的各种统计平均量,
对于任意格点 Gx ∈,任意组态 S∈a,记
}(min)( }\{}\{ xGxUxxGxU HmH x ahag h
DD
==,(9,7)
那么
),(}{ baxp ∑
∈
=
Gy
mH
mH
xyGyU
xxGxU
e
e
])([
])([
}{\
}{\
ab
ab
=
≥=≥ eGeGGe x
xGxGUU HH
||
1
||
1
||
):)()(max( }{\}{\
d
gbgb
,
其中 xd 为 能量函数 UH 在 格点 x 上的振幅,?为 UH 的 振幅,其 确切 定义为
235
|)()(|max
}{\}{\
bad ba UUx HH
xGxG
= =,xGx d∈?=? max,(9,8)
由此 推出
D≤?≤ epGPC
xx 1),(min||1)( }{,}{ baba,(9,9)
从而 得到
||
}{}{ )1()()()( 1
G
xx ePCPCPC M
≤= L ||1 Ge≤,(9,10)
定理 9,１２ （ 遍历 收敛定理） 以上定义的 转移矩阵 P的 Dobrushin 收缩系数满足
||1)( GePC≤ (?为
UH 的 振幅 ),
记以 P为转移矩阵的 Markov 链为 nx,则它以 Gibbs 分布 p 为可逆分布,而且有
P pTnn 1 →? ∞→,
证明 显见 ),()(),()( abbpbaap JJ pp =,由 P Sp ∈= baba,)),(( 的定义立得
),()(),()( abbpbaap pp =,因此 p 为可逆分布,又由 1)( <PC,用 Dobrushin 定理可知极限成立,
[ 注 ] G - 格点的排序 ),,{ ||1 Gxx L 并不是实质的,我们可以用一个 G 上取值全正的预选概率分布 )(xg,( 即 1)(,0)( => ∑
∈Gx
xgxg 例如均匀分布 ) 的 随机数作 随机选取来代替,这时转移矩阵 P应作相应的改动,代替 ),(}{ ba
kx
p 的是,以 )(xg 的概率对只在 x 坐标不同的组态间进行转移,即令
P (?= ~P ||) G,
~P
Sp ∈= baba,
~
)),((,
其中
==
)(0
)(),()(),( }\{}\{}{~
其它情形若 xGxGxpxgp bababa
,
这时 我们仍有
P pTnn 1 →? ∞→,
１,５ Gibbs 分布的模拟退火
本段 讨论 求 得达到 组态空间 上 能量函数最小值的组态 的概率方法,就是用能量函数的
Gibbs 分布作模拟退火的方法,考虑依赖于温度 nT 的 Gibbs 分布,以它们作为 不变分布构造各步 为 非时齐 Gibbs 转移 矩阵，由此 确定的 一个 非时齐 M arkov 链 }{ nx,我们可以用它作模拟退火,使得当 n 充分大以后，此非时齐的 Markov 链 nx 以 大概率落在 能量函数最小的组态
236
所组成的 集合 中,
更 具体地说,对于 如下的 依赖于温度 nT 的 邻位势 Gibbs 分布
)(1)( 1
)(
a
ap Un
HT
n
n e
Z
=,∑
∈
=
S
HT
n
U
neZ
a
a )(1
,
相应地定义
=?
)(0
),()(1),(
)(1
,
)(
\
其它外坐标相同在 JeZp JSJUn
HT
Jn
n
J
baaba
ab
,(9,11)
∑?=
J
JSJU
n
HT
Jn eZ
b
ab
a
)(1
,
\)(
,
再令
P ()( }{
D
=nx Snxp ∈= baba,)( }{ )),((,P =)(n ∏
∈Gx
P )( }{nx,(9,12)
则 P )(n 是一个转移矩阵,我们 将 以它为 非时齐转移矩阵 列 所确定的非时齐 M arkov 链记为
}{ nx,
定理９,１３ （ Gibbs 分布 的模拟退火的收敛性定理 ）
如果温度 nT 满足 如下的 Geman - Geman 条件,对 00 ↓< nT,存在 0n 使 0nn ≥ 后有
n
GT
n ln
||?≥,(9,13)
其中 D 是能量函数 )(aUH 的振幅,
)(min)(max aaD aaD HH SUS ∈∈?=,
那么在 任意 初 始 分布 0m 下,以 P )(n 为第 n 步转移矩阵列的非时齐的 Markov 链 }{ nx,在时刻 n 的绝对分布 nm 有极限
)(∞→ pm
n,
其中 )(∞p 为集合 min})(:{ =∈ aa UHS 上的离散均匀分 布,从而 还 有
1))(min)((lim == ∈∞→ ax a USnUn HHP,
证明 我们验证 Dobrushin - Isaacson - Madsen 定理的条件成立,条件 (A.1) 的验证与第
8 章 3,2 段的模拟退火定理中的证明完全一样,下面验证 Dobrushin 条件 (A,2) ’,这里定理９,１２中 的不等式变为
237
≤ nT
G
n ePC
||
)( 1)(,
其中?为能量函数 UH 的振幅,于是对于 )()1()( njj PPP L+,我们有
≤+ )( )()1()( njj PPPC L ≤+ )()()( )()1()( njj PCPCPC L )1(
||||
nT
GG
jk
e
=
∏,
由 无穷乘积的理论,当 ∞→n 时,上面右方趋于 0 的充要条件为级数 ∑
n
T
G
ne
||
是发散的,
但是由 Gemam - Gemam 条件,我们有
∑
n
T
G
ne
||
∑?≥
n
ne ln ∑ ∞==
n n
1,
所以 Dobrushin 条件 (A,2)’成立,
［注 1 ］ Gibbs 场的模拟退火 在 图象的清污 ( 图象的滤波 )中的应用
设图象为组态 a,而在传输过程这受到了污染,使得接受方观测到的图象为组态 b,需要用数学 方法进行 滤波,以 便 恢复图象的清晰度,我们把在观测为组态 b 的条件下,状态 a 的分布,记为 )|( baP
(就是后验分布 Posteriori distribution,即 Bayes 分布 ),一般它也可用写成 Gibbs 分布的形式,
)|(
)(
1)|( ba
bba
He
ZP
=,(9,14)
其中 )|( baH 称为 条件能量函数,典型的例子如
)|( baH ∑ ∑
∈
+?=
)(
2
2
2 ))()((
2
1))()((
yx x
xxyx basaaq,
其中第一项强调了图形的光滑性,第二项为保真度,滤 波的目的 就是要找组态
^a
,使条件能量函数达到最小,
)|(min)|( ^ baba a HH =,
这个最佳组态
^a
可以 通过 模拟退火方法求得,
又例如,在边缘检测中,b 就是离散化点上的象素,而 a 就是,边随机链,的一个样值,与滤波问题的不同的是,此时 的 后验图形只是先验 图形的一部分,即边的位置,在复杂的图形处理技术中,例如,X 断层摄影术的处理,情况就远 为复杂,适当的利用分割技术可以得到较好的恢复效果,
[ 注 2 ] 图象的处理有许多数学方法,Gibb s 分布方法只是其中的一种,另一种更为常用的方法是,把图象的一些特征参数作为观察量 ( 这同时也达到了图象信息压缩的目的,用它们来估计 Ma rkov 链 的状态,
就是 重建图象,这就是隐 Markov 方法 (HMM),一般地说,如果特征参数取得合适,HMM 方法是十分有效的,
例如,纹理图象分析方法 (Texture Analysis) 提取特征参数,其它的方法 还 有,下节将介绍的随机函数迭
238
代系统方法,除此以 外,还有小波分析等非概率方法,
[ 注 3 ] 当格点集 G 为无限集时,组态的数目就不再是可数的了,此时就应归入连续状态的 Markov
场的范畴,在统计物理中应用的模型,就是这种模型,有时自旋空间 S 也可以是连续的,例如,圆周 或其它集合,这时就可能存在多个不变分布,称为有相变,因为每个分布代表了统计物理中质点的一个宏观
“相,,统计物理中也正是对这种相变的存在性更感兴趣,
[ 注 4 ] 若每次发展的修改 的 格点集合 J 不是单点集,且各次的修改集可 以 不同,那么,此时的 Markov
链未必有遍历定理,这 时并 不能保证此 Markov 链 nx 的时间平均收敛到不变分布,因此,即使 n 很大,也不能认为 nx 的样本近似为不变 分布的样本,
2 时间离散状态连续的 Markov 链
2,1 概率空间再访
设 }:{ 基本事件?=? ww,F 是样本空间? 的某些子集组成的集合类 (事件体 ),满足
对于任意 ∈nAA,F )1( ≥n,恒有 IU
∞
=
∞
=
∈?=
11
,,
n
n
n
n AAAA W
D
F,
( 这样的 事件体 F 代表了可以通过理论推算,或实际测量能够得到概率的全体 随机 事件 ),在概率理论中 F 称为? 上的一个?s 代数,(?,F ) 称为 概率空间,于是概率 P 就是定义在
s 代数 F 上的一个非负函数,
定义９,１４ （ Borel 集与 Borel 函数）
设 B 是 dR 的某些子集组成的一个?s 代数，且满足以下条件,
(1) dR 中的 任意开集 ∈G B,
(2) 任意一个有 dR 的某些子集组成的?s 代数 G,只要它含有一切开集,则就有
B? G,
那么 B 称为 dR 上 Borel 集合类,其中的集合称为 Borel 集,即 Borel 集合类是以开集为元素的最小?s 代数,( 以前在本书中的所谓,常见的,集合,均指 Borel 集 ),1R 上的 Borel
集称 为一维 Borel 集,
再设 M 为 dR 上的一个满足下述条件的函数类,
(1) 对极限封闭,若 ∈nf M 且 )()( xfxf n →,则有 ∈f M,
(2) 任意连续函数 ∈f M,
(3) 若 dR 上一个函数类 L,只要对极限封闭,而且包含全体连续函数,则有 M? L,
239
那么 M 称为 dR 上 Borel 函数类,其中的函数称为 dR 上的 Borel 函数,( 以前在本书中的所谓,常见的,函数,均指 Borel 函数 ),
再则,随机变量的概念是依赖于 事件体 (?s 代数 ) F 的,严格地应称为 F - 随机变量,
定义９,１５ 样本空间? 上的实值函数 x 称为 F - 可测的,如果对于任意实数 a,恒有 ∈≤ })(:{ awxw F,于是,x 为 F - 随机变量与 x 为 F - 可测的是一样的,而随机变量的分布则依赖于概率 P 的給法,
分布函数为 )( xF 的 随机变量 x 的函数 )(xg 的 数学期望 是如下的 Stieltjes 积分
∫= )()()( xdFxgEg x ( def= ∫ )()( dxFxg )
( 括号中的 )(dxF 就是 )(xdF,前一种记号便于推广到多变量,或条件分布函数 的情形 ),
而
)|)(()|()|( ],( hxhxh xIExPxF?∞
=≤= = hhx =?∞ = yx yIE )]|)(([ ],(
称为 x 对于 h 的 条件分布函数,类似地,我们有
∫= )|()(]|)([ hhx dyFyggE,
( 注意,右方的积分是直观地理解的,如果想要一个清晰的证明,在数学上还必须要作一系列,后台操作,),
2,2 时间离散状态连续的 Markov 链
定义９．１６ 实值随机序列 }0:{ ≥nnx 称 为（时间离散状态连续的） Markov 链,如果对于任意 0≥n,任意实数 10,,,,?nxxxy L,有
====≤+ ),,,|( 00111 xxxyP nnnn xxxx L )|( 1 xyP nn =≤+ xx,(9,15)
(这正说明条件分布函数有 Markov 性）．在 ),,,( 01 xxx Lnn+ 有联合密度时,(9.15)等价于 1+nx
的条件分布密度的如下的等式
====
+
),,,|( 0011
1
xxxyp nnn
n
xxxx L )|(
1
xyp n
n
=
+
xx,(9.15)’
与第５章中类似地，这里的 Markov 性可以有如下的多种等价性描述,
Markov 性的等价性质 1
对于任意实 Borel 集 10,,,?nLLL L,以及随机事件 },,{ 1100 ∈∈= nnA LxLx L,恒有
)(),( 11 xPxAP nnnn =∈==∈ ++ xLxxLx,(9,16)
240
当存在联合密度时,(9.16)式也有其相应 的条件分布密度形式,
Markov 性的等价性质 2
对于随机序列 }{ nx 在大于时刻 n 所确定的随机事件 },,{ 11 knknnnB ++++ ∈∈= LxLx L
及随机事件 },,{ 1100 ∈∈= nnA LxLx L,有
)(),( xBPxABP nn === xx (9,17)
或
)()|()( xAPxBPBAP nn === xx,(9,17)’
（ 直观 证明 我们只在一切条件分布密度都存在情形证明 (9.17)式,先看 2=k 情形,由条件概率的性质,推广了的全概率公式与 Markov 性的定义,我们有
(),( PxABP n ==x ),|,2211 Axnnnnn =∈∈ ++++ xLxLx
（直观地） ),,|( 1122 AxP nnnnn =∈∈= ++++ xLxLx ),|( 11 AxP nnn =∈ ++ xLx
dyyfAxyP nnnnn
n
)(),,|( 1122
1 ++++Λ
==Λ∈= ∫
+
xxx )|( 11 xP nnn =Λ∈ ++ xx
dyyfyP nnnn
n
)()|( 1122
1 ++++Λ
=Λ∈= ∫
+
xx )|( 11 xP nnn =Λ∈ ++ xx,
其中 1+nf 是随机变量 1+nx 的分布密度,由于上式右方与随机事件 A 无关,可见等式左方应该等于
)( xBP n =x,对于一般的 k,只需做归纳法,）,
Markov 性的等 价性质 3
在已知“现在”的条件下，“将来”与“过去” 是条件独立的 (此即 (9,17)’)
Markov 性的等价性质 4
对随机序列 }{ nx 及 0,1 ≥≥? nm 及,任意 Borel 集 L 及实数 10,,,?nxxx L,都有
====∈+ ),,,|( 0011 xxxP nnnmn xxxLx L )|( xP nmn =∈+ xLx,(9,18)
（ 直观 证明 设 一切条件分 布密度都存在,m=1 时即为 Markov 性的定义,对 m 作归纳法,假定 m
时 (9.18)正确,今证 m+1 时它也正确,
),,,|( 00111 xxxP nnnmn ===∈++ xxxLx L
dyyfxxxyP nnnnnmn )(),,,|,( 1001111 ++++ ====Λ∈= ∫ xxxxx L
),,,,|( 001111 xxxyP nnnnmn ====Λ∈=+++∫ xxxxx L
dyxxxyf nnnn ),,,|( 00111 ===+ xxx L
利用归纳法假设及 Markov 性的定义,上式简化为
)|( 11 yP nmn =∈= +++∫ xLx dyxyf nn )|(1 =+ x,
241
注意上式 右方与 },,{ 0011 xxnn == xx L 无关,可知 (9.18)式对 1+m 也成立）,
Markov 性的等价性质５
对于任意有界 Borel 函数 g,即任意实数 10,,,?nxxx L,均有
))((),,,)(( 1100 xgExxgE nnnnn ==== +?+ xxxxxx L１,(9,19)
(在 )()( xIxg Λ= 时,(9,19)就成为 (9,15),然而，从 (9,15)推导 (9,19) 需要测度论的基本知识,本书从略．这里我们粗略地指出其朴素思想,就是?+?= ggg,其中 }0{ >+ = ggIg,而 ±g 们可用有限个示性函数的线性组合来近似 ),
注 等价性质 1,2,3,5 诸款都有相应于等价性质 4 的形式，读者可自己列出它们,
Markov 性的等价性质 6 （最一般的形式）
对于任意实 Borel 集 Λ,及只由随机序列 }{ mx 在时刻 n 及其后的信息所决定的随机变量 nh,恒有
)(),,( 11 xPxxP nnnnnn =Λ∈===Λ∈ xhxxh L (9,20)
或
)()),,,( 0011 xExxxE nnnnnn ===== xhxxxh L,(9,21)
注 这个等价条件的内涵十分丰富,其等价性的证明在测度论中非常典型，本书从略,
与离散状态的 Markov 链的转移矩阵类似,在连续状态 情形,我们有概率转移核,
2,3 概率转移核
定义９．１７ 记
=),(),( Lxp mn )(),|( nmxP nm ≥=∈ xLx,(9,22)
这是一个条件分布 (上式右方就是 )|)(( xIE nm =xxL ),
∈==
)(0
)(1)(),(),(
L
LdL D
L x
xxxp nn,
特别地,),()1,( Lxp nn + 称为时刻 n (到时刻 1+n )的 概率 转移核,而 ),(),( Lxp knn + 就称为从
n 出发的 k 步概率转移核,
特别地,若存在 )(),,(),( mnyxp mn <,使
∫= LL),(),( xp mn dyyxp mn ),(),(,(9,23)
那么 ),(.),( yxp mn 就是,在 xn =x 的条件下,mx 的条件分布密度,称为 条件转移密度,
242
（我们在本书中如无特别声明,则恒假定 条件转移密度 存在．）它满足以下性质,
(P.1) )(1),(0 ),( mnyxp mn <≤≤,∫ 1),(),( =dyyxp mn,
(P.2) ( Chapman - Kolmogorov 方程 ) 对于 nm ≥≥? 有
=),(),( yxp ln ),(),( zxp mn∫ dzyzp lm ),(),(,(9,24)
有时也简写成 乘积核 的形式,即记成
=),( lnp?),( mnp ),( lmp,
( (P.1)是显然的,(P,2)的直观含义十分清楚,它就是条件概率的全概率公式 ),
(P,2)的证明 由 Markov 性及条件期望的性质
)|)|)((( nmlIEE xxxL )|),|)((( nnmlIEE xxxxL= )|)(( nlIE xxL=,
从而
),(),( zxp mn∫∫
Λ
dzdyyzp lm ),(),( )|),(( ),( xdyypE nmlm == ∫
Λ
xx
)|)|)((( xIEE nml == xxxL )|)(( xIE nl == xxL dyyxp ln ),(),(∫
Λ
=,
（注 当不存在条件分布密度时,Chapman - Kolmogorov 方程应写为
=),(),( Lxp ln ),(),( dyxp mn∫ ),(),( Lyp lm ),(9,25)
其中积分是 Stieljies 积分）
2,4 时齐的连续状态 Markov 链
定义９．１８ 连续状态的 Markov 链称 为 时齐的,如果其概率转移核 )()1,( L，xp nn + 与
n 无关 (即 1 步转移概率与出发时刻 n 无关,我们把它改记为 ),( Lxp,
在时齐且有条件转移密度的情形,我们简记
),(),( )1,0( yxpyxp
def
= ),()1,( yxp nn +=,),;,0(),()( ymxpyxp m?= ),(),( yxp mnn +=
这时 Chapman - Kolmogorov 方程就简化为
∫=+ dzyzpzxpyxp nmmn ),(),(),( )()()(,(9,26)
在本书中除非特别声明，所考虑的连续状态 Markov 链总是 时 齐的,
(注 在不存在 条件 转移密度时,由 (9,24) 得到
=+ ),(),( Lxp mnn?+ )1,( nnp（?++ )2,1( nnp ),(),( )(),1( LL xpxp mmnmn =? +?+ ）L
它也不依赖出发时刻 n,它是 Markov 链 }0,{ ≥nnx 的 m 步转移核,其中
∫= ),(),(),()2( dyxpypxp LL,
∫∫ == ),(),(),(),(),( )1()1()( dyxpypdyxpypxp mmm LLL,)
在本 书中除非特别声明，所考虑的连续状态 Markov 链总是 时 齐的,
243
连续状态的 Markov 链的条件转移密度 (或不存在密度时的转移核 ),是刻画此链的基本统计特征量,一个 Markov 链的统计性质可由其初始状态 0x （或其初始分布
)()( 00 xPx == xp ）以及其条件转移密度 (或不存在密度时的转移核 )完全确定,这就是下面的定理
定理 9,１９ （ Markov链的有限维分布密度） 若连续状态的时齐 Markov链 }0,{ ≥nnx
的初始状态为 x=0x,其条件转移密度为 ),( yxp,则
(1) 在条件 x=0x 下 ),,( 1 nxx L 的条件联合分布密度为
),(),(),( 1211 nn xxpxxpxxp?L,
如果 0x 具有分布密度 )( xj,那么,),,( 0 nxx L 的联合分布密度为
)( xj ),(),(),( 1211 nn xxpxxpxxp?L,
(2) 对于有界 n 元 Borel 函数 g 有
== )|),,([ 01 xgE n xxx L
nnnn dxdxdxxxgxxpxxpxxpn LLLL 2111211 ),,(),(),(),(12?ΛΛΛ∫ ∫ ∫,
（证明提示 (1),对于 n 用归纳法,(2) 对于 )()( 11 nn xgxgg L= 可用归纳法证明,而对于一般的 g,则需要用测度论中的典 型逼近方法,这里从略,(其粗略思想为,用示性函数的线性组合近似
}0{ >
+ =
ggIg,,再用形如 )()( 11 nxIxI nΛΛ L 的函数的线形组合,去近似一般类型的示性函数
),,( 1 nxxI LΛ,其中 L 是任意 n 维 Borel 集 ),
注 在转移核不存在密度且 0x 没有分布密度 的情形,),,( 0 nxx L 的联合分布为
∫∫
ΛΛ
=Λ∈Λ∈
n
nnP LL
0
),,( 00 xx )(0 dxF ),(),(),( 1211 nn dxxpdxxpdxxp?L,
其中 )(0 xF 是 0x 的分布函数,),( Λxp 是 Markov 链的转移核,而
== )|),,([ 01 xgE n xxx L
),,(),(),(),( 11211
12 nnn
xxgdxxpdxxpdxxp
n
LLL?ΛΛΛ∫ ∫ ∫,
定理９,２０ （连续状态的 Markov 链的绝对概率 ）
设 )|( xyf n 是在条件 x=0x 下,随机变量 nx 的条件分布密度．那么
244
dzzxpzyfxyf nmmn ),()|()|( )( ）（∫=+,(9,27)
证明 应用全概率公式,
2,5 例
例９．２１ ( 非随机的迭代映射 ) 设随机 变量 序列 nx 满足迭代关系
xf nn ==+ 01 ),( xxx,
这时 nx 没有分布密度,事实上它是 一个常数,即由 x 只能转移到 )(xf,于是,对于任意实
Borel 集 L 可以直接计算得到
),,,|( 00111 xxxP nnnn ===∈+ xxxLx L
),,,|)(( 0011 xxxxfP nnn ===∈= xxxL L })({ L∈= xfI,
也就是 其转移核为
)(),( )(1 xIxp f LL?=,
其中 )(1 L?f 是 L 在 f 下的原象集,即 })(:{)(1 LL ∈=? yfyf,
例９．２２ ( 带随机噪声的迭代映射所生成的连续状态的 Markov 链 )
设 }{ nw 为独立序列，且 nw 的分布 密度 为 )( yf n,),( yxg 是 2R 上 的 Borel 函数．又设随机变量序列 nx 满足迭代关系
),(1 nnn wg xx =+,
其中 0x 与 }{ nw 独立,此时 nw,),( nwxg 都与 },,,{ 10 nxxx L 独立．故而
),,,|( 00111 xxxyP nnnn ===≤+ xxxx L
),,,|),(( 0011 xxxywxgP nnnn ===≤= xxx L
∫ ≤=≤= dzzfywxgP nyzxgn )()),(( ),(,
上式右方与 01,,xxn L? 无关,从而等于 )|( 1 xyP nn =≤+ xx,可见 nx 是连续状态的 Markov
链,而且,当且只当 )( yf n 不依赖 n 时,这个连续状态的 Markov 链才是时齐的,
2,6 寻找 dR 上可微函数 )(xf 的最小值的位置的模拟退火算法
先考虑 1=d 情形,设 nw 为独立同分布的 )1,0(N 随机变量序列，又设 函数 0)( ≥xf
245
且二阶光滑，再设 )(xf 取最小值的点全体组成集合 0S,我们考虑如下的随机迭代模型
001 ),0()()(' xhfwhf nnnnnn =>++=+ VVeVVV,
它是例９．２ ２的特殊情形，而且在 xn =V 的条件下,转移到 yn =+1V 的 条件转移 密度为
=+ ),()1,( yxp nn
2]
)(
)('[
2
1
)(2
1 xf xhfxy
n
ne
xf
e
pe
,
即它是正态分布 ))(),('( 22 xfxhfxN ne+,其中 x 作为参数,于是 nV 是时间离散的非时齐的状态连续的 Markov 链,这里的非时齐性是由于“小参数” ne 依赖于时刻 n 而引起的,
对于 1>d 情形,则相应的随机迭代模型为
001 ),0()()( xhfwfh nnnnnn =>+?+=+ VVeVVV,(9,28)
其中 nV,nw 都是随机向量,而 f 为多元函数 ( 注意如果 0Sn ∈V,则有 0)( =? nf V,此时随机迭代模型就变成
001 ),0()( xhfw nnnnn =>+=+ VVeVV,
即只有随机小扰动项 )( nnn fw Ve 起推动 nV 运动的作用,由此可知一旦当此 Markov 链的轨道 nV 到达了 f 的某个稳定点以后,再到它的给定的某个小邻域以外的可能性就较小，而且当此稳定点是 f 的最小值点时，这种可 能性就相对地更小,我们再注意,当 此 Markov 链的轨道 nV 在未到达 f 的某个稳定点之前，在（９．２８）中的项 )( nfh V? 起的主要作用将使
nV 往 f 的稳定点发展．综合此两个朴素思想的启示,我们可以用此 Markov 链的发展来“近似”地找到 f 的最小值的位置，这正是连续 状态空间变量的模拟退火思想的要点,当然用这种方法还是有失败的可能的．然而，这样做最终失败的可能性较小,如果我们能忍受小概率失败的失误而放弃,而只求以较大的概率能达到整体最小值附近,那么这是一种可行的妥协方法 ( 这种加入随机小扰动项 )( nnn fw Ve 的方法，能避免用确定性方法 ( 例如梯度方法 )
求最小值位置时,陷入局部极小处的失误 ),把 ne 看成温度后,（９．２８）提供了连续状态空间的一种模拟退火方法,即利用增加小扰动项的随机迭代模 型作模拟退火,如果适当地选取退火的温度 ne 趋于 0 的速度（例如,如 Geman - Gemam 那样选取 1)(ln?= nn ge,而 g 又足够大）,那么,在理论上可以证明,对于任意 0>d 有
1))(( 0 →?∈ ∞→nn SBP dV,
246
其中 )( 0SBd 指 0S 中所有点的δ邻域的并 ( 这 种极限关系称为 nV 按 概率收敛到集合 0S,与有限状态空间的情形看起来不同,在 有限状态空间 情形 是收敛到最小值集合,而现在时收敛到最小值集合的一个邻域,如果我们认为有限集中的点都 是孤立的,那么这两者是一致的 ),
也就是说,对随机差分方程作迭代,经过足够多次以后,就以近似于 1 的概率,保证 nV 会到达 f 取最小值的集合的较近的范围,这就 是用随机迭代模型作模拟退火的基本思想，这里温度 ne 依赖 n 的阶 1)(ln?= nn ge,对于概率收敛性的成立是至关重要的,
注意求函数 )(xg 零点的问题,可以化为求 2)(xg 的最小值问题,因而也可以应用随机迭代模型的模拟退火方法,因此这个数学模型在广泛的领域中有应用前景,
[ 注 ] 模拟退火方法中的人工随机干扰 nw,是 独 立同分布随机 序列,称为 白噪声,模拟退火的进程 一般较 慢,对此 人们采取了 多种加速收敛措施,其中常用的有 TINA ( Time
Invariant Noise Annealin g,用时间不变的噪 声 作 模拟退火 ) 算法,它的基本思想是,采用,有色噪声,作为逃逸出局部极值 的,驱动力,,即在目标函数值大的时候,使干扰大些,
而当函数值接近最小值时,就把干扰取得小些,特别是在已知最小值的情形 ( 例如求解非线性联立方程组时,只要把目标函数取为各方程 的 平方和,那么目标函数的最小值就是 0,
而 达到最小值的变量就是该方程组的解 ),例如,在已知最小值的条件下,对 )(xf 求整体最小值 的位置,也 可以考虑如下随机差分方程,
nnnnn wfff ))(()( *1?+?+=+ xxxx,(9,29)
其中 nw 是独立同分布的随机变量列,*f 是目标函数 )(xf 的最小值,于是 (９,２９）给出了用有色噪声 nn wff ))(( *?x 干扰的 TINA 算法,
我们还可以将 TINA 算法修改为适合于未知最小值的全局优化问题,例如可以把干扰噪声取为 nn wncfg ])ln()))(([ 1?+?+? gx,其中 g 是一个有界增函数,c 是 一个常数,
它取得越接近最小值 *f 越好,而在实际计算中,人们可以根据对 f 的 函数值的认识,不断地调小参数 c,以求加速收敛,
2,7 Dobrushin 不等式,指数遍历性与收敛性
在第 5 章中，我们曾对离散状态的 Markov 链,讨论了当 n →∞时它的极限性 态,在本节中，我们将利用 Dobrushin 不等式，对于连续状态的 Markov 链得到相应的结果．这类结果要求假定转移密度存在，并满足较强的条件，这时不仅可以得到 n 步转移密度当 n →∞ 时的极 限趋于唯一的不变概率分布密度，而且保证误差具有指数型的下降速度 ( 称 之 为指数遍历性 ),这一结论比之于一般 Markov 链的遍历定理要强得多,
命题 9,２３ ( Dobrushin 不等式 )
对于 Borel 函数 f,定义 dxxff |)(||||| ∫=,假定时间离散的连续状态的 Markov 链
247
存在转移密度 ),( yxp,那么对于任意分布密度 )(),( 21 xx rr 有
||||)(||)(),()(),(|| 2121 rrrr?≤ ∫∫ PCdxxxpdxxxp,(9,30)
其中
dzzypzxpPC yx |),(),(|sup21)(,?= ∫ （９．３１）
称为转移密度 ),( yxp 的 Dobrushin 数,】
例９,２４ 设 )(xf 是一个有界连续函数．那么对于满足 )()( yfxf ≠ 的点 ),( yx,
2
)()(),( yfxfyxh +=L 是方程
2
))(( 2xfz
e
2
))(( 2yfz
e
=
的唯一解,且 ),( yxh 也是有界连续函数,当 )()( yfxf < 时我们有
∫ ≤? dzzypzxp |),(),(| ),( yxh∞?∫ p21 dze
xfz
2
))(( 2 ∞
∫+ ),( yxh p21 dze
yfz
2
))(( 2
)(),( xfyxh?
∞?∫
＝ p21 due
u
2
2
∞
∫
+
)(),( yfyxh p2
1 due u22?
2m
∞?∫
≤ p21 due
u
2
2
∞∫+
1m p2
1 due u22? ))(1()(
12 mm Φ?+Φ=,
其中
}|:)(),({|inf,1 yxxfyxhm yx ≠?=
D
,}|:)(),({|sup,2 yxxfyxhm yx ≠?=
D
,
∫
∞?
=
a
a)(F p21 due
u
2
2
,
显见 ∞<<<∞? 21 mm,于是
yxPC,sup2
1)( =
p2
1 2 ))(( 2| xfze∫ dze yfz |2 ))(( 2
2
1≤ 1))](1()([
12 <Φ?+Φ mm,
例９．２５ 设连续状态的 Markov 链的转移密度 ),( yxp 满足条件
0)(),( ≥≥ ygyxp,
那么此转移密度的 Dobrushin 数满足
248
1)(1)( ≤?≤ ∫ dyygPC,
证明 仿照第５章离散状态的 Markov 链的 Dobrushin 不等式的证明可知
=?= ∫ dzzypzxpPC yx |),(),(|sup21)(,
= dzzypzxpzypzxpyx )),(),((sup ),(),(,?∫ ≥ dzzgzxpzypzxpyx ))(),((sup ),(),(,?≤ ∫ ≥
≤ Sup dzzgdzzgzxpyx )(1))(),((,∫∫?=?,
定义９．２６ 具有转移密度 ),( yxp 的连续状态的时齐 Markov 链 nx 称为 指数 L 1遍历 的 ( 或 强指数遍历 的 ),如果存在正整数 0n,分布密度 )( y）（ ∞r 及常数 0,10 ><< cr,只要 0nn ≥,就对任意 x 满足
nn crdyyyxp <? ∞∫ |)(),(| )()( r,
其中 ),()( yxp n 是此 Markov 链的 n 步转移密度,
下面的概念只需要转移核，而不需要假定存在转移密度,
定义９．２７ 连续状态的时齐马氏链 nx 称为 指数遍历 的，如果存在正整数 0n,分布函数 )()( yF ∞ 及常数 0,10 ><< cr,只要 0nn ≥,就对任意 yx,满足
n
n cryFyxP <?≤
∞ |)())((| )(x,
这里 改 用 )(xnx 记初值为 x=0x 的 Markov 链 nx,分布函数 )()( yF ∞ 不依赖于初值 x,称为
nx 的 ( 遍历 ) 极限分布,也称为 nx 的 平稳分布,以
)(∞F 为初分布，就能得到一个连续状态的平稳过程,显见,指数 1L 遍历一定指数遍历,
定理９．２８ （指数 L 1遍历定理）
具有转移密度的连续状态时齐 Markov 链 nx,如果存在 10 ≥n,使其 0n 步转移密度
),()( 0 yxp n 的 Dobrushin 数 1)( )( 0 <nPC,那么此 Markov 链是指数 L 1遍历的,
( 证明 与第 5 章中离散状态情形 类似,从略 ),
例９．２９ 例９．２５中 的 Markov 链,只要 存在 10 ≥n,使其 0n 步转移密度满足
),()( 0 yxp n )( yg≥,且 ∫ > 0)( dyyg,
则它 指数 L1遍历
249
例９．３０ 设 f 是有界函数，而 nw 是独立同分布随机变量序列，且具有严 格正的连续的分布密度 )( xp,那么由如下定义的迭代列 nx 所决定的连续状态 Markov 链是指数 L 1
遍历的
nnn wf +=+ )(1 xx,(9,32)
证明 不妨设 bxfa ≤≤ )(,此 Markov 链的转移密度为
min))((),( ≥?= xfypyxp bza ≤≤ )( zyp? ＞ 0,
由例９,２５和定理９．２８可知此 Markov 链是指数 L1遍历的,
定理９,３１ （指数遍历定理）
设 nw 是独立同分布随机变量序列,),( vuf 是连续函数，且关于 u 单调不减,又
)),(()(,10 nnn wxfxx xxx == +,(9,33)
如果存在 z 0,0n 及 0>d 使
))((sup 0
0
zxP nx ≤x 及 ))((inf 0
0
zxP nx ≥x 都 d≥,
则连续状态的 Markov 链 )(xnx 是指数遍历的,
( 证略 ),
3 随机的迭代函数系统 ( 随机 IFS,随机 Iterative Function System)
3,1 局部相似性的基本想 法
如第２节所示，图象的局部相似性是它的基本特点,最简单的图象是只有黑白两个灰度的轮廓图,由于这种图的轮廓线在局部总有一点相似性 ( 所谓自相似性 ),人们就希望用函数迭代映射的,不动点,,来近似描绘轮廓线图,并进一步把一般的黑白灰度图和彩色图用随机迭代函数映射的不变分布来近似,这样就可以达到在传输过程中压缩图象信息,而在接收后再重建合成人工图象的目的,( 显然这是非常朴素的想法,例如,压缩图象信息的压缩比如果能达到 20,则效果就很不错了,但是在实践中,还必须考虑到与原有的设备能相容这个方 法,即改进压缩比的方法只希望在原有设备上附加一个系统,而不希望完全改建设备,因此再理想的压缩方法也必须与实际情况相协调 ),
图像的纹理分析 (Texture Analysis)是一种重要的图像处理方法,在目标识别等方面有重要应用。图象的纹理是一种结构,它包括：自然纹理表面的灰度，颜色，条纹,以及它们的模拟,要模拟或分析图像纹理，选择合适的模型是关键的一步,纹理模型必须能表述图像元素的空间分布特征,均匀的纹理的最基本的分布特性是它的空间相似性，就是在图像上的不同地方的纹理有着几何上相似的象素组 合,在数学上，可以用小波分析,自相似随机过程,
随机迭代函数系统,Markov 随机场等来描述纹理,由于 Markov 随机场和 Gibbs 随机场是等价的，因而就可以用 Gibbs 随机场模型 (用另一种看法视为随机图,而把图像作为随机图的
250
一个样本 ) 描述纹理．进而可以使用各种方法估计 Gibbs 分布的参数，在此基础上就可以对纹理进行分析和处理,下面介绍的轮廓图的方法,是一种与 Gibbs 分布不同的,近似描绘 无灰度的黑白图的方法，即迭代函数系统方法,
3,2 轮廓图全体组成的距离空间
定义９,３２ （轮廓图的图象窗口中的距离）
只有轮廓而无灰度的黑白图,即轮廓图,其黑色部分就可看成为平面上某个窗口的一个闭集,我们记此窗口为 X,它是 2R 中的有界闭集,例如,]1,0[]1,0[ ×=X,X 中的两个点 ),(),,( 2121 yyyxxx == 之间可定义距离 ),0[,∞→× XXd,例如普通的 Euclid 距离 222211 )()(),( yxyxyxd?+?= ( 也可以是其它的距离,只要满足如下的 3 条距离公理,
(D.1) 0),( ≥yxd,而且等号成立的充要条件为 yx = ;
(D.2) ( 对称性 ) ),(),( xydyxd = ;
(D.3) ( 三角形不等式 ) ),(),(),( zydyxdzxd +≤,
例如,|)||,max(|),( 2211 yxyxyxd=,或者 ||||),( 2211 yxyxyxd?+?= ),
轮廓图的图象空间与图像之间的距离
如上所述,一幅轮廓图就是窗口的一个闭子集,从数学上描述两幅轮廓图的近似程度的最朴素的想法是,在两幅图之间引入距离 ( 也可以引入准距离,相对熵,,这将在第 10 章的隐 Markov 模型与第 15 章的算法中叙述 ),为此,我们先介绍一些关于距离空间的基本知识,由前面的论述,我们知道,窗口 X 可配以距离 d,我们把 ),( dX 作为蓝本,引入距离空间的概念如下,
定义９．３３ 任意一个集合 Y,如果定义了满足 3 条距离公理的映射 ),0[,∞→× YYd,
则我们称 ),( dY 为 距离空间 (metric space),距离空间中有收敛概念,对于 Yxn ∈,如果有 0),( →xxd n,则称 nx 收敛到 Yx ∈,记为 xxn →,
距离空间中的序列 )(},{ Yxx nn ∈ 称为 Cauchy 列,如果对于任给 0>e,存在 N,只要 Nmn >,,就有 e<),( mn xxd ( 显见这个概念是欧氏空间中的 Cauchy 序列的推广,需要注意的是,在距离空间中并没有加法运算,除非它还具有线性空间结构 ),
距离空间 ),( dY 称为 完备的,如果它的任意 Cauchy 列 }{ nx 都有极限,即存在 Yx ∈,
使 0),( →xxd n,( 例如,dR,dR 中带边的超立方体,dR 中带边的有界图形,都是完
251
备的距离空 间 ),
对于距离空间 ),( dY,Y 的子集 F 称为 闭集,如果 F 中的任意收敛序列的极限仍在 F
之中,闭集的余集称为 开集,Y 的子集 K 称为 紧集,如果 K 中的任意序列都有收敛 子序列,
而且此收敛子列的极限仍在 K 之中,
显见紧集必是闭集,在欧氏空间中,紧集与有界闭集是一样的,在上述的窗口空间中,
由于窗口本身是有界的,因而闭集与紧集是一样的,但是,在一般距离空间闭集就未必是紧集,
定义９,３４ （作为距离空间的轮廓图的图象空间）
设 ),( dX 为窗口距离空间,在此窗口空间上的一个轮廓图,就是窗口空间的一个闭集,
记
}:{)( XFFXH?=? 闭集,
那么,)( XH 为所有的轮廓图的图象组成的集合,称为 图象空间,对于任意两个轮廓图
)(,XHBA ∈,可以定义如下三种不同的距离 21,,rrr,
),(sup),(sup),( AydBxdBA ByAx ∈∈
+=r,
其中
),(inf),( yxdBxd By∈= ;
},:{inf),(1 ddd dr ABBABA=?,
其中
}),(:{ dd <∈?= yxdAxyA 使 ;
)],([sup)],([sup),(2 BydBxdBA ByAx ∈∈
∨=r,(9,34)
其中
),max( baba =∨,
不难验证 21,,rrr 都是 )( XH 上的距离,而且它们都彼此等价,即收敛性是一样的 ( 虽然距离的大小并不一样 ),
命题９．２５ 距离空间 )),(( rXH 是完备的,
证明 设 nA 为 )),(( rXH 的 Cauchy 列,则易见对于任意 nn Ax ∈,nx 是 ),( dX 的 Cauchy 列,
故有极限 x,令
})(:{ xAxxA nn →∈?=?,
用与普通数学分析类似的方法,可以证明 A是闭集,即 )( XHA ∈,而且 0),( →AAnr,这就证明了
252
)),(( rXH 的完备性,
完备距离空间的压缩映射与不动点
定义９．２６ 任意距离空间 ),( dX 自身的映射 W,如果存在 1<c,使得对于任意
Xyx ∈,,恒有
),())(),(( yxcdyWxWd ≤,(9,35)
那么,W 称为 压缩映射,c 称为 压缩系数,
例９．２７ (仿射压缩映射 ) 设 W 为如下映射
+
→?
=
2
1
2
1
2221
1211
2
1
a
a
x
x
aa
aa
x
xx W,
假定系数
2221
1211
aa
aa,
2
1
a
a 满足：
1
1,
0
1,
1
0,
0
0 的 W 像都在 X 中,且
=c
2221
1211
aa
aa T
aa
aa
2221
1211 的最大特征值 1<,那么,可以证明 这样的仿射映射是压缩映射,
定理９．２８ ( 不动点定理 ) 对于任意完备距离空间 ),( dX 上的压缩映射 W 有
(1) 存在唯一的不动点 aaWa =)(,,
(2) (迭代近似 ) 记,2 WWW oo = )),(()(2 xWWxW =o,1 nn WWW oo o=+
则对于任意 Xx ∈,恒有 axW n →)(o,
(3) (近似不动点 ) 如果 y 满足 e<)),(( yyWd,则 cayd?≤ 1),( e,
证明 先证明不动点的存在性,令 =nx )(xW no,于是 ),(),( 11 xxdcxxd nnn ≤+,从而
∑
=
+++ →?≤≤
p
k
n
knknnpn xxdc
cxxdxxd
1
11 0),(1),(),(,
所以 }{ nx 是 Cauchy 列,由完备性,它有极限,记为 a,再则
0),()),(())(),(()),(( →++≤ axdxxWdxWaWdaaWd nnnn,
故 aaW =)(,再证明不动点只有一个,假 设 另有不动点 b,则 有
),())(),((),( bacdbWaWdbad ≤=,
253
由 1<c 立得 0),( =bad,即 ba =,这就证明了 (1)和 (2),最后
)),(())(),(())(,(),( 1 ayWdyWyWdyWydayd nnn oooL +++≤?
)),((1 ayWdcc nn oL ++++≤? eee c?≤ 1 e,】
从窗口空间 ),( dX 的压缩映射导出图象空间 )),(( 1rXH 的压缩映射
定理２．２９ ( Collage 拼贴定理 )
若窗口空间 ),( dX 上有有限个压缩系数为 ic 的压缩映射 iW )1( Ni ≤≤,定义
)),(( 1rXH 上的映射
^W
,其象集为所有的 )( NiWi ≤ 的象集之并,
)( XHA ∈?,U
N
i
i AWAW
1
^
)()(
=
=,
那么,
^W
是 )),(( rXH 上的压缩映射,
^W
称为由 },,{ 1 NWW L 导出的 迭代函数系统,它的唯一不动点 0A 称为 迭代 函数 系统
^W
的 吸引子,而且对于任意 )( XHF ∈ 有误差估计 ( 称为拼贴定理 ) ( 设 ic 们中的最大数为 c )
))(,(1 1),( ^101 FWFcAF rr?≤,(9,3 6 )
证明 由 )( XH 上的等价距离 1r 的定义得到
))(),(( ^^1 BWAWr U U
N
i
N
i
ii BWAW
1 1
1 ))(),((
= =
= r
))(),((max 1 BWAW iii r≤ ),(1 BAcr≤,】
迭代函数系统的吸引子的＂自相似＂特性
由 },,{ 1 NWW L 导出的迭代函数系统
^W
的吸引子 0A 满足,
UN
j
ji AAWAWAW
1
00
^
00 )()()(
=
==?,
这说明 吸引子 0A 的 iW 象是 0A 的一部分,再用 kW 们继续地迭置作用下去,则它相当于对
0A 不断地作细分,在国内外工程文献中,称以上性质为,迭代函数系统
^W
的吸引子 0A 关于 },,{ 1 NWW L 有自相似性,事实上，当 kW 们都是仿射压缩映射时，这种说法较为合理,
而在一般情形是一种近似的描述,在实际运作中，迭代函数系统中的映 射常取成仿射压缩
254
映射，这就充分反映了轮廓图的自相似性,
用迭代函数系统压缩图象信息的思路为：通过选取 N 及 },,{ 1 NWW L,用迭代函数系统
^W
的吸引子 ( 在距离 1r 下 ) 去近似标的图象．由此可以达到压缩信息的目的，因为这时轮廓图的信息完全可由 },,{ 1 NWW L 描述,下面的定理说明上述要求不但能够做到,而且可以把 iW 们选取为仿射映射,
定理９,３０ ( Hutchinson 定理 )
任给的轮廓图 )( XHF ∈ 及容许误差 0>e,存在 1>m,及压缩的仿射映射组
},,{ 1 mWW L,使
^W
的吸引子 0A 是 F 的容许 近似,即
er ≤),( 01 AF,
证明 因为 F 是紧集,所以可以用有限个半径为 8e 的圆 mBB,,1 L 覆盖 F,U
m
k
kBF
1=
,我们不妨假定这些圆的中心都在 F 的 8e 邻域中,于是 4
1
e
FBF
m
k
k
=
U,我们可以找压缩系数 21<ic
的仿射映射 iW,使 ii BFW?)(,又由于 iB 的直径小于 4e,故
2)()(
e
FWBFW iii,
从而
4
^ )( eFFW?
,2
^ )( eFWF?
,
因此，由距离 1r 的定义得到
2))(,(
^
1
er ≤FWF
再用 Collage 定理推出
er ≤),( 01 AF,
[注 ] Hutchimson 定理出现于 1981 年,它的主要含义是,可以把轮廓图近似地作为迭代函数系统的吸引子,吸引子具有,自相似性,,也就是工程界广泛认为的,分形性质,,
轮廓图是窗口 X 的一个闭子集，可以用迭代函数系统的吸引子来建模．而窗口上的灰度图的相对灰度可以看成口上的一个概率分布,为此考虑用随机迭代函数系统生成的
Markov 链的不变分布来建模．这就是下一段的理论,
在窗口 ]1,0[]1,0[ ×=X,相对灰度图就相应地考虑为 概率 空间,( X B ))( X 上的一个概率分布,其中 B )( X 为 X 的 Borel 子集全体组成的类,
255
3,3 灰度图与随机迭代函数系统
本段考虑最简单的随机映射 及由它 生成的 Markov 链，并用以建模．这种随机映射就是,对于 有限个映射，以不同的概率从中取一个．其确切的定义为
定义９．３１ 窗口 ]1,0[]1,0[ ×=X 自身的 随机压缩映射 W的定义为,存在 X 自身的压缩映射 )()1(,,mWW L,使
(P W )1(,) 1)( =++== mii pppW L,
直观地,我们可以 将 它 记为 如下的分布 列的形式 （与离散分布的不同处是其取值为映射，而不是数）
W
)()1(
1
,,
,,~
m
m
WW
pp
L
L,(9,37)
我们 称
)()1(
1
,,
,,
m
m
WW
pp
L
L 为 随机映射 W的分布,
命题９,３２ （独立同分布的随机映射生成的 Markov 链）
设 LL,,,1 nWW 为独立同分布的随机映射序列,设 ０x 是随机变量．定义
xxx,),( 011 LW= =+1n W )(1 nn x+,(9,38)
则 }{ nx 是时齐的 Markov 链,称为 随机迭代函数系统,其概率转移核 (不存在条件转移密度 )
为
),( Lxp ∑
=
=
m
i
Wi xIp i
1
)( )(1)( L,(9,39)
证明
),,,|( 00111 xxxP nnnn ===∈+ xxxLx L
(P= W ),,,|)( 00111 xxxx nnnn ===∈+ xxxL L
(P= W ))(1 L∈+ xn ∑
=
∈=
m
i
i
i xWPp
1
)( ))(( L ∑
=
=
m
i
Wi xIp i
1
)( )(1)( L,
推论９．３３
∑
=
+ ==
m
i
i
inn xWfpxfE
1
)(
1 ))((]|)([ xx,】
下面介绍独立同分布的随机压缩映射生成的 Markov 链不变分布的存在与唯一性,
定理９．３４ ( 遍历性定理 )
随机压缩映射生成的 Markov链 nx 存在唯一的不变分布,设它是某个随机变量 V 的分
256
布函数 )xF（V,那么对于此 Markov 链的转移核 ),( Lxp,分 布函数 VF 有 以下的不变性质,
=)(xFV )(]),(,( ydFxyp V∫?∞,
其中 ],(],(],(),,( 2121 xxxxxx?∞×?∞=?∞=?,VF 称为 不变分布函数,而且,对于任意有界 Borel 函数 f,不管什么初值 0x,恒有遍历定理，即有如下的按概率收敛式
)()()()()(1
1
∞→=?→?∑ ∫
=
nxdFxfEffn
n
i
p
k VVx,
特别地,如果把 nkk ≤≤1}{x 中落入 ],( x?∞ 的频数记为 ],( nn?∞,则有以下的按概率收敛式,
)()(],( ∞→?→∞ nxFnn px V,
这就得到了不变分布 函数 的相合估计如下 (如果将一幅灰度图用一个随机迭代函数系统 建模,
那么,利用 此相合估计 可以在实际计算中,通过随机迭代函数系统重建灰度图 )
=)(
^
xFV nn x],(?∞,(9,40)
证明 令 x== 00 xh,定义
hhhhh,),(),( 12011 Lnn WW == + =+1n W )(1 nh,
那么 }{ nh 与 }{ nx 同分布,由于 }{ nh 是某个概率空间,(? F,)P 上的随机变量序列,对固定某个基本事件 (样本 )?∈w,存在 mii n ≤≤,,1 L１ 使 )()( )()( 1 xWW niin oLo=wh (其中运算 o 表示映射的复合 ),于是对于 此固定的基本事件 w,压缩映射 )()1(,,mWW L 中最大的压缩系数 )1(<c 有
))(),(( whwh npnd + )(( )()( 1 xWWd nii oLo= ))((,)()( 1 xWW pnii +oLo
0))(,( )()( 1 →≤ ++ xWWxdc pnn iin oLo,
这说明 ))({ whn 是 X 中的 Cauchy 列,从而它收敛于某个 (依赖于初值 x 的 ) )(wV x,即存在随机变量
xV (理论上可以严格证明 ∈≤ }{ axV F ),使
1)( = →? ∞→ xnnP Vh,
再则,由于 )(whn 也依赖于初值 x,我们 将 它改记为 )(wh xn,又 因为它是压缩映射的复合,所以仍然是压缩映射,故 有
257
0),())(),(( →≤ yxdcd nynxn whwh (关于 w 是 一致 的 ),
由此得到
))(),(( wVwV yxd ))(),(( whwV xnxd≤ ++ ))(),(( whwh ynxnd ))(),(( wVwh yynd 0→,
这说明极限随机变量 xV 与初值 x 无关,所以我们可以 将 它改记为 V,显见,nh 的分布函数
n
Fh 按分布 收敛 到 V 的分布函数 VF ( 即在 )( xFV 的连续点 x 上有 )()( xFxF
n Vh
→ ),
由 nx 与 nh 同分布得到,nx 的分布函数
n
Fx 也 按分布 收敛到 V 的分布函数 VF,可以证明,这个分布 VF 就是转移核 ),( Lxp 的唯一的不变分布 函数,
(注 此 事实可以证明如下,对于任意有界连续函数 f,有 (f x ()1 fn =+ W )(1 nn x+ )
(Ef W ）＝)(1 nn x+ ∑
=
m
i
n
i
i WEfp
1
)( ))(( x ∑
=
=
m
i
n
i
i WfpE
1
)( ))](([ x
∫∑
=
= )())((
1
)( xdFxWfp
n
m
i
i
i x ∫∑?
=
= ))(()( 1)(
1
xWdFxfp i
m
i
i nx,
但是按定义有 (f x ()1 fn =+ W )(1 nn x+ ),对此 两边取期望便得到
)()(
1
xdFxf
n +∫ x ∫ ∑
=
= ))(()( 1)(
1
xWFpdxf i
m
i
i nx,
令 ∞→n,利用 按分布 收敛的性质 推出
)()( xdFxf V∫ ∫ ∑?
=
= ))(()( 1)(
1
xWFpdxf i
m
i
i V,
由 f 的任意性,再利用（ 9.39） 式便 得到
)( xFV ∑
=
=
m
i
i
i xWFp
1
1)( ))((
V )()(
1
],()( 1)( ydFyIp
m
i
xWi i V∑ ∫
=
∞?= )(]),(,( ydFxyp V∫?∞=,
(把它 写成测度的形式就是
=)(LVF )(),( ydFyp VL∫ ),
这就证明了极限随机变量 V 的分布 VF 是转移核 ),( Lxp 的不变分布 函数,再则,如果此转移核另有不变分布 函数 F,将 以它为分布 函数 的随机变量 记为 J,取 Jx =0,那么,由 迭代定义的随机变量序列
}{ nx 就是一个平稳 序列,与 上面的证明 类似地 仍有 1)( =→ VxnP,于是作为极限分布的 F 应该等于
258
VF,此外,还是由于 1)( =→ VxnP 及 平稳列的性质,还 可以 进一步 证明有不依赖初值的遍历定理 ),】
(对于具有一些抽象积分知识的读者,我们 还 可以用距离空间的方法来证明不变分布 函数 的 存在 性如下,
令 M 为窗口空间 ),( dX 上全体概率测度的集合,即 M 中的 每 一个元素 m 都是 X 上的 Borel 集合类
B )( X 上定义的一个概率,在 M 上可以定义如下的 Hutchinson 距离 r,
∫ ∫?= ∈? |)()()()(|sup),( xdxfxdxfLf nmnmr,
其中 )},(),,(|)()(:|{ XyxyxdyfxffL ∈?≤?=?,可以证明 ( M,r ) 是完备距离空间,对于压缩系数为 ic 的压缩映射 )(iW 们组成的随机映射
)()1(
1
,,
,,
m
m
WW
pp
L
L
,我们定义它在 M 上导出的 一个 映射为 (记为 1?W )
∑
=
=→
m
i
i
i WpW
1
1)(1 )( mmm,
往证 1?W 是压缩的,
= ),( 11 nmr WW ∫ ∫∈? |))(()())(()(|sup 11 xWdxfxWdxfLf nm
∫ ∫∑?=
=
∈ |)())(()())((|sup
)()(
1
xdxWfxdxWfp ii
m
i
iLf nm
∫ ∫∑?≤ ∈
=
|)())((1)())((1|sup )()(
1
xdxWfcxdxWfcpc i
i
i
i
Lf
m
i
ii nm
∑
=
≤≤
m
i
ii cpc
1
),(),( nmrnmr,
这里我们用了,对于 Lf ∈,有 LWfc i
i
∈)(1 o,而 c 为 ic 们中的最大者,一旦证明了压缩性,那么 1?W 作为完备距离空间 ( M,r ) 的迭代映射的不动点就存在唯一,而且可以用迭代逼近,这个不动点就是我们所要的不变分布 函数 ),
遍历定理中关于压缩映射的条件可以减弱为 如下的 平均压缩性,这就是 Barnsley 与
Elton 在 1988 年得出 的 平均压缩条件下的遍历定理
定理９．３５ 若随机映射
)()1(
1
,,
,,
m
m
WW
pp
L
L 满足平均压缩条件,存在 1<c,使对于任
259
意 Xyx ∈,有
∏
=
≤
m
i
p
ii yxcdyWxWd
i
1
),())](),(([,(9,41)
那么,对于由独立同分布的如上的随机映射生成的随机的迭代函数系统，遍历定理的所有结论仍然成立,
(证 明 略 ),
(以上的不动点 (不变分布函数 ),作为图 像有某些自相似性,Barnsley 等人进一步希望用某种不动点来描述自相似的图象,为此提出了递归迭代函数系统 (Recursive Iterative Function System),他用一个 Markov
转移矩阵 P 来代替随机迭代函数系统中的随机向量 ),,( 1 mpp L,即记以 P 为转移矩阵的时齐的 Markov
链为 }{ nh,再 定义
)()(1 nn nW xx h=+,
此时 }{ nx 不再是 Markov 链,但是 },{ nn hx 却是 Markov 链,同样可以得到它的不变分布 函数,
另一方面,我们还可以在窗口上的全体概率测度的集合上,定义一个压缩映射得到不动点,而且它与
上面定义的是一样的,进一步还可以考虑有多个递归迭代函数的系统,得到这个系统的不动点 ),
4 统计中的 Bayes 方法与 图象的处理,分割 与重建
4．１ Bayes 统计要义
Bayes 方法的基本想法是,把概率函数 (包括分布密度和离散的分布概率 ) 中的 未知参数 q 当作随机变量 (或随机向量 ),
在 Bayes 方法中不再区分参数与随机变量．所以未知参数的分布的确认是最重要的．在抽取样本以前，就只能根据先验知识设置未知参数的分布，称为先验分布，或验前分布．在抽取样本以后，根据对于抽取到的样本的概率规律的了解，就可用 Beyes 公式把参数先验分布改进为后验分布，或验后分布，也称 Bayes 分布,
１．先验分布与后验分布 (Bayes 分布 )
设随机变量 x （或随机向量）的分布（或概率函数）为 ),( qxp （ q 可以是向量）．用
Bayes 的观点，把 q 看成随机变量后,),( qxp 就是：在 q 已知的条件下,x 的条件分布，
即
),()|(| qqqx xpxp =,
假定 q 的分布密度（或概率函数）为 )(Jf,称为 q 的 先验分布 (a priori distribution),,
于是 ),( qx 的联合分布为 ),()( qqf xp,在 J 取固定值的条件下,假若 nxx,,1 L 为采自
),( qxp 的独立随机样本,在 nxx,,1 L 已知的条件下,q 的条件分布记为
260
),,|( 1 nxxJj L,由 Bayes 公式可知它是 q 与 ( nxx,,1 L )的联合密度 ),()(
1
qxqf i
n
i
p∏
=
，除以 ( nxx,,1 L )的边缘密度，即
),,|( 1 nxxJj L
∫ ∏
∏
=
=
n
i
i
n
i
i
dp
p
1
1
),()(
),()(
JJxJf
JxJf
,(9,42)
称为 q 的 后验分布 ( a posteriori distribution),或 Bayes 分布,
2,已知损失函数是的的 Bayes 估计
假定对参数 J 的一个估计
^J
带来的损失可由一个 损失函数
^),( JJL
来度量．我们用
BayesE 表示在后验分布下取的数学期望,对于 J 的任意一个估计量
^J
,它在后验分布下的平均损失 ),(
^
JJLEBayes 称为 后验风险,使后验风险达到最小的估计,称为 Bayes 估计,记为
Bayes
^J
。即它满足
)],([min)],([ ^^ ^ JJJJ
J
LELE BayesBayesBayes =,( 9,43 )
3,Bayes 方法评述
Bayes 方法的长处是设置了一个先验分布，它可以发挥对未知参数的已有统计知识的认知．然而其难点也正在于如何选取合适的先验分布,这可以说是这种由从另一个角度认识 统计带来的一对矛盾,在实际应用中正是要巧妙地处理好这一对矛盾,
再则，由于在 Bayes 方法中，参数与随机变量是平等的，所以对于随机变量的估计也常用 Bayes 估计．特别 是 在数据不全的情形（缺失数据）的估计中，需要用测量数据来估计缺失数据，或者是在滤波问题中，或隐 Markov 模型中，需要用量测数据估计状态 （参见第
10 章与第 15 章）,
4．先验分布的参考取法
第一种取法出现在无任何先验信息的情形．具体地说
（ 1）在 J 取值于 ),( ∞?∞ 时，认为它是均匀的（这称为 Bayes 假设）,即先验密度为常数 (如果 q 可能取值区域面积为无限,则它不是真正的分布,而是所谓的,广义,分布，但是此时 (9,
43)仍然是一个分布，仍称为后验分布或 Bayes 分布 ),那么,Bayes 估计就变成为普通的最大似然估计,
（ 2）在 J 取值于 ),0[ ∞ 时，常用广义分布 J1 。
（ 3）在 J 取值于 ]1,0[ 时，常用广义分布 )1( 1 JJ? 。
第二种取法是用 最大 熵原则：取所谓的最大熵分布（参见第 10 章），它是在已知的限制条件下，信息掌握得最少的分布，说一句白话,就是“最吃亏的”分布,在对于参数完
261
全没有限制的情形下,广义的 最大熵分布恰恰就是均匀分布,所以第二种取法也是第一种取法的延伸。又若参数是一个概率向量（或者要估计的是一个取值于,概率向量组成的概率空间,的随机向量,这里
,n 维概率向量组成的概率空间” }1,1,,0:),,{( 111 =++≤≤=? nnn pppppp LLL ），
那么这时的最大熵分布是 Dirichlet 分布。这就是在生物信息论中常用 Dirichlet 分布作为先验分布的原因,
经验知识对于粗估先验分布常常很有用,例如,若未知的是英文字母或汉字,人们可以根据字母或汉字出现的经验频率得到先验分布的统计估计,
第三种取法是采用 Jeffreys 原则：取 先验密度为 )(det)( qqf IC=,其中 C 是常数，
)(qI 是信息矩阵，它的分量 )),,((),),(ln()( 1
2
k
ji
ij
xpEI qqq
qq
qq L=
=,
第四种取法是统计中的所谓共轭分布方法,如果参数的先验分布与它关于分布 ),( qxp
的后验分布属于同一种类型，那么此种类型的分布就称为是 ),( qxp 的 共轭分布,人们认为用共轭分布作为先验分布有其合理性，因为这时先验分布与后验分布属于同一种类型,
常见分布的共轭分布有,Beta 分布是二项分布 的共轭分布,也是几何分布 的 共轭分布 ;
Gamma 分布是 Poisson 分布的共轭分布，也是指数分布的共轭分布；正态分布是已知方差的正态分布的共轭分布； 双侧截尾正态分布是已知方差的双侧截尾正态分布的共轭分布 ; 逆
Gauss 分布 是已知 均值 的正态分布 ( 包括 双侧截尾 情形 )的共轭分布 ; Perato 分布
（ )0,(),()( ),[)1( >?= ∞+? axIxaxp a gg ga ）是均匀分布的共轭分布 等,
在对参数没有任何先验知识时，在上面四种取法中，到底用那个更为合适？一般地说，
并没有确切的答案．这是 Bayes 方法的不足 之 处．实际上这时只能随 问题的不同作具体的考虑,
5,Bayes 估计再访
通过 Bayes 分布来估计参数或状态，有 3 种常见的方法,
（１）对给定的损失函数，求使 后验风险最小的估计,
（２）用后验分布的数学期望,即 后验期望估计,
（３）用后验分布达到最大值的点来估计，称为 最大后验估计,
与一般数理统计类似地还可以用后验分布作区间估计,
6,Bayes 策略
定义９,３６ 在统计抉择中，通常会涉及样本，参数与行动．参数 q 是被考虑成随机变量的,将 可能采取的行动全体组成的集合 记为 A．假定在对象的状态分布的参数为 q 时，
如果采取行动 )( Aa ∈,就会导致 损失 ),( aL q,而采取的动作 a 是依赖于被抽取到的样本
nxx,,1 L 的,我们将 它记为 ),,((),( 1 naa xxxx L==,这个“函数”（从样本到动作的对应） )(?a 就称为一个 策略 。 在样本为 x 时，策略 )(xa 带来的损失为 ))(,( xq aL,它 是样本的函数,
262
定义９,３７ 在 q 固定的时的平均损失为
xdxpxpxaLaELaR n ),(),())(,())(,())(,( 1 qqqxqq L∫==,(9,44)
称为 策略 )(?a 关于 q 的风险函数 。它对于先验分布 )(qf 的平均
=?))((ar ∫? qqfq daR )())(,(
(9,45)
称为 策略 )(?a 对于先验分布 f 的 Bayes 风险 。使先验分布 f 的 Bayes 风险达到最小的策略，
称为 对于先验分布 f 的 Bayes 策略 。
定义９,３８ 动作 )(xa 带来的损失 ))(,( xq aL 关于 q 的后验分布的期望，则称为 策略
)(xa 的后验风险 （就是 动作 )(xa 的后验风险），其表达式为
qxqxqxJ dpaLaLEBayes )|())(,())(,( ∫=,(9,46)
其中 )|( xqp 为参数 q 的 后验密度,即后验分布的密度,
这里需要分清楚,策略的风险函数是对状态变量作的平均,即对于一切可能的样本作平均，所以它定义在整个策略类上，而策略对先验分布的 Bayes 风险，是策略的风险函数关于先验分布的平均．而后验风险是，在给定样本下，策略对参数的后验分布作的平均．从数学处理上,在求后者时样本是固定的，所以要简单得多．这就是使用后验分布的好处，当然，后验分布的计算也常常并不简单,
需要注意，,Bayes”的冠名，通常用来修饰后验性质，例如,Bayes 分布就是后验分布．但是，在 Bayes 风险的定义 中，却被用来冠名先验性质,这似乎在冠名上有些混淆，
但是因为有下面的定理，所以并无实质上的不妥,
定理９．４０ (Bayes 策略基本定理 )
一个策略 )(*?a 是 Bayes 策略的充要条件是,对于任意样本 x,)(* xa 都使后验风险达到最小,
（此定理的证明需要较深的知识,从略）,
这个定理说明了,Bayes 策略在任意样本 x 上的取值，就是 x 的最小后验风险动作．于是最小 Bayes 风险在本质上仍反映了后验性质,
由此定理可见，只要算出了后验风险,再求得 Bayes 策略就并不困难。
Bayes 估计可以纳入 Bayes 策略的框架．再者，假设检验也可以纳入 Bayes 策略的框架，
因而 Bayes 统计方法的重点就归结为 Bayes 分布的计算。
在实际问题中,Bayes 分布的计算，常常要通过随机模拟得到,特别是 Bayes 分布中的分母的计算，常常涉及巨大的工作量，较为有效的数值计算方法，就是 Markov 链 Monte Carlo
方法（ MCMC）,
4,２ Bayes 方法在图像中的应用与 观测量不是状态变量时的参数估计
图象处理包括,分割 ( Segmentation),图象压缩,图象重建,图象分类 (如利用相似性对纹理图分类 ),图像识别,聚类等等,
263
图象在传输过程 先 需要进行压缩,然后由接收方 再将 压缩图象信息恢复,所以,图象处理主要为,图象压缩与图象恢复,而图象分割常用在,例如,自然场景的分割,在图象恢复,
图像分割的各种方法中,Bayes 方法是很重要的一个,它归结为在把接收到的信息作为已知的 条件下,对描述图象的一些参数的重估,
Bayes 聚类方法,就是根据观测 (从传输渠道得到 ),在轮廓图 (或灰度图 )全体组成的集合中,选取一个在此观测的 条件 下,Bayes 风险最小的轮廓图 (或灰度图 ),作为恢复的图象,
这里损失函数
^),( JJL
的取法是十分重要的,它必须取得使估计的结果较为稳健 (robust),
即更多地与分布的定性行为相关,而对具体分布的表示式的敏感性尽量地小,这里不能一味简单地搬用模式识别中的现成结果,而必须考虑到,标的图象,的来 源,结合它所在学科的特点,具体情形具体分析,
图像的先验分布 )(ap (Priori distribution),其中 a 为组态,一般可以写成 Gibbs 分布的形式 (只要它处处为正 ),
)()( aap He?=,
其中
)(ln)( apa?=H,
称为 能量函数 或 验前质量指标,
图像的信息特征中,一般还有代表边缘信息的 ‘边能量,( 对应于,边过程,),各部分的能量强度 的不同常表现为某些加权参数 jq 的大小,例如,图象恢复的模型可以取
∑
=
=
3
1
}|()|(
j
jjHH baqba,
这是一个参数建模问题,就是要确定参数,即要确定 H 为能量函数的 Gibbs 分布的参数 jq
们,在观测 )()1(,,mbb L （接收到的信息）条件下,作最大似然估计或 Bayes 估计,此时的联合分布 ( 即似然函数 ) 为
∏
∑=?
∑
∑
=
=
=m
k H
H
j
kj
j
k
j
e
eL
1 )|(
)|(
321 3
1
)(
3
1
)(
),,(
g
bgq
baq
qqq,
这里既有待估参数 q,又有未知组态 a,要同时估计它们,一种有效的方法是将在第 15章中介绍的 EM 算法,即用 EM 方法作迭代近似,对于格点数太多的情形,Geman 兄弟提出了一种简化的修正算法,可以大大地减少计算量,其基本思想为,一旦只要得到目标指向局部极大的
)(^ n
q 后,就按 ),|(
)(^
bq
n
P? 取样
)(^ n
a,然后按下面的方法更新
)1(^ +n
q,
任取
)1(^ +n
q 满足 (P
)(^ n
a ),|
)1(^
bq
+n
),|(
)(^)(^
bqa
nn
P>,
264
在实际估计中,对于参数的先验知识的利用是至关重要的,例如,参数所在的位置或区域,
Gibbs 场的 Bayes 分布纹理分析时,主要是根据纹理的特征,选取合 理的 Gibbs 分布,即选取能量函数,常用的能量函数有 Graffigne 的?Φ 模型,
∑?
+
=
相邻),(
2
2||
1
1)(
yx yx
xyH
d
aaqa,
在模型选取好以后,纹理图的聚类就化为：通过观测到的一系列标准样品,进行学习训练以确定 (估计 ) 参数 q,再用 Bayes 聚类或其它聚类法来决定此纹理图形属于那一类,
图像识别的融合原则
在图象识别时,有时可以利用不同的识别器的融合 (fusion),这时可以采取的原则,
例如有,多数原则,权威原则,民主原则 (以相似度 作 加权 )等,
习题 9
1,设 nw 独立同分布,且具有正连续密度 )( xr,又 bxga ≤≤< )(0,
)1(mod)(1 nnn wf +=+ xx,nnn wf +=+ )(1 hh,那么连续状态 Markov 链 nn hx,都是指数 1L 遍历的,
2,设 nw 独立同分布,且具有正连续密度 )( xr,假定 0>yf,且存在 0)( ≥yg,使在 x 固定时
),( yxf 对 y 的反函数,记为 ),(1 yxf?,满足 )(),(),(( 1 ygy yxfyxf≥?r,而 Markov 链 nx 满
足 ),(1 nnn wf xx =+,那么 nx 是指数 1L 遍历的 Markov 链,
3,设 nw 为独立同分布列，具有分布密度 r,证明
110011,,,),,(+?+ ==+= kknknnn xxwf xxxxx LL
是 k 阶 Markov 链,意即
),,,|( 0011 xxxAP nnnkn ===∈+ xxxx L )|( xAP nkn =∈= + xx,
龚光鲁，钱敏平著 应用随机过程教程及其在 算法与智能计算中的应用
清华大学出版社,2003
第 9 章 以图象信息为背景的随机场,迭代 Markov 系统
随机场的一个直观背景,是以图像为基本事件的随机元素,另一方面,迭代 Markov 系统是取值为压缩映射 随机映射 的迭代,其目的是得到其 不变元素,以用来描述一个图像,
1 有限格点上的 Markov 随机场与图象
1,1 有限格点上的 Markov 随机场
定义９,１ （ 相邻系统与子团 ）
设 G 为 平面有限格点 集合,}1,1:),{( 21 NjNijiG ≤≤≤≤=,对于 G 中的任意两点 yx ≠,可以定义它们相邻,或者不相邻,这样定义了 S 中的一种相邻概念以后,就得到了 G 上的一个 相邻系统,注意,两点是否相邻完全可以是抽象的,并不需要在直观上比邻
( 有时为了看起来有对称性,也可以认为 1N 与 1 恒同,2N 也与 1 恒同,即认为 G 为环面格点,)
在本书中,我们把 yx,相邻记为 yx?~,记 与 x 的相邻 的 点的 全体 为 )(x?,G 的子集
C 称为 G 的 子团,如果 C 中的任意两个元素都相邻,含 x 的子团 记为
U )()( xxC?=,
一个相邻系统由 G 中点对的相邻概念所决定,因此 也由全体子团所决定,所以人们常用? 来代表一个相邻系统,或用全体子团 }{C 来代表一个相邻系统,
例 ９．２ (1) 若 f=? )(x,则只有单点集是 G 的 子团,
(2) 若 }{\)( xGx =?,则 G 的所有子集都是 子团,,
(3) 一阶相邻 ( 最紧相邻,简称 紧邻 ),? 的相邻点 o可以 如下图
o
oo
o
此时 G 的全部 子团 形式 为,
f,?,,
|
(4) 二阶相邻,? 的相邻点 o可以 如下图
230
ooo
oo
ooo
此时 G 的全部 子团 形式 为,
f,?,,
|,
\ 及其转动,
\| 及其转动,
×
||,
高 阶相邻等 可以 如 下 图 所 示
6
54345
42124
631136
42124
54345
6
其中数字代表分阶 为 相邻的阶,S 的一个相邻系统也就是以 G 的点为顶点的一个双向图,
定义９,３ （ 随机场 ）
以 平面矩形格点集 G 中的格点 为参数,取值于 }1,,0{?LL ( 也可以是 实直线 R,或一个离散的集合 ) 的随机变量族,称为 取值于 }1,,0{?LL 的 G － 随机场,简称为 随机场,
按照 统计物理中的术语,我们称 }1,,0{?LL 为 自旋空间,
于是,随机场的一个样本,就是在格点集 G 的每个格点上,坐着,}1,,0{?LL 中某 一个 数,随机场的样本的一个直观 背景是图像．黑白图象是连续变化的一个灰度图,而 在实际图象处理中的 黑白灰度都是经过离散化采样的,设黑白灰度可分成 L 个 (一般 L=256)不同的等级,}1,,0{?LL 中的一个数值可以解释为黑白的灰度,最粗放的是黑白轮廓图，在每个采样点上只有,黑,,“白,两个象素之一 (即 1=L 的情形 ),一般 G 为 256× 256 个采样点 (即 25621 == NN ),那么,一幅黑白灰度图就可以看成在 G 的每一个格点上的一个灰度 ( 象素,Pixel),即给每一个格点赋予 1,,0?LL 中的一个值,所以,G － 随机场的样本,
即做一次试验可能得到的结果,可以解释为 G 上的一个黑白灰度图,
在这个意义上,随机场就也可以解释为 随机图,或者说,它是在图像空间取值的一个随机元素 ( 简称为 随机元 ),这是随机场的第 2 种解释,从相对灰度 ( 每个格点上的灰度除以总灰度 ( 灰度的和 )) 的角度也可以说,随机场的分布是格点集 G 上的一个随机的概率分布,
定义９,４ （组态） 从数学的角度抽象,随机场的一个样本就是 GLS }1,,0{?=? L 的一个元素,记为 Lba,,称为一个 组态 ( configuration ),即 ):( Gxx ∈= aa,
231
}1,,0{?∈ Lx La 称为组态 a 的 x 坐标,再记 ))(:( GJxxJ?∈= aa,称为组态 a 的 J
坐标,
所以,随机场也可以看成 是一个 随机组态,或者看成 为在 组态空间 GLS }1,,0{?= L
上取值的随机元,这是随机场的第 3 种解释,从分布的角度 也可以说,随机场 的 分布 是 组 态空间上的一个概率分布,
由此可见,一个灰度图就可以 近似地 视为 一个组态,
定义９,５ （ Markov 随机场 ）
随机场 }:{ Gxx ∈x 称为 Markov 场,如果对于任意格点 Gyx ∈,及任意组态
GLS },,1{ L=∈a,恒有
),|( xyP yyxx ≠== 一切axax )~,|( xyP yyxx
=== 一切axax,(9,1)
Markov 场 的含义是：在 任意格点 x 的余格点位置 上取 已知 值的条件下，随机场在格点 x 处取值的概率规律只与格点 x 的? 相邻点有关,
由定义 可以 直接 证明 下述定理
定理９．６ 随机场 }:{ Gxx ∈x 是 Markov 场 的 等价 条件为,对 任意格点集
GA? 及任意组态 GLS }1,,0{?=∈ La,恒有
),|,( AyAxP yyxx?=∈= axax ))(,|,( AyAxP yyxx?∈=∈== axax,
其中 }{)( 中点的邻点全体AA?=?,
1,２ 相邻系统 的 Gibbs 分布与 Gibbs 随机场 ( 邻位势 Gibbs 场 )
定义９．７ （ 格点集 G 的 子集与组态的相互作用 ）
对于 给定 的 一个子 格点集 )( GA?,组态上的一个 与之相关的 函数 )(aAU
( GL }1,,0{?∈ La ) 称为 A 对于组态 a 的相互作用,如果 它 满足,对于任意 两个组态
GL }1,,0{,?∈ Lba,只要 )()( bxax
xx = ( Ax ∈? ),就有 )()( ba AA UU = ( )(aAU
代表组态 a 在 A处的,坐标,对 a 的作用,所以 要求对在 A 处,坐标,相同 的 任意 组态的相互作用 一 样 ),
定义９．８ （ 位势函数 ） 在补充定义 0=fU 后,相互作用
}:{ GAU A?
的 全体 称为 位势,对于给定的相邻系统?,如果对于任意 A，只要它不是 子团，就有
0=AU,那么 }:{ 子团的为GCUC 就 称为 邻位势,如果 相邻正好就是直观上的 紧邻,则这样的 的 邻位势 简称 紧邻位势,
232
定义９,９ （ 能量函数 ） 对于 邻位势,组态的函数
∑
=
C
CU UH
子团
(9,2)
称为 邻位势的能量函数,
定义９,１０ （ 邻位势 的 Gibbs 分布与 Gibbs 随机场 ）
在组态空间 GLS }1,,0{?= L 上定义的概率分布 p ∈= aap (),( S ),称为 邻位势 ( 在能量函数下 ) 的 Gibbs 分布 ( 这里是离散分布 ),其中
)(1)( aap UHe
Z
=,
UH 为 邻位势 的能量函数,而 Z 为一个归一化常数,称为 Gibbs 分布 的 配分常数,
∑
∈
=
S
HUeZ
a
a )(,(9,3)
以 Gibbs 分布 为分布 的随机场 称为 Gibbs 场,即
)(),( apax =∈?= GxP xx,
就是说,Gibbs 分布 p是随机场在所有的格点处取值的联合分布,
在 第 6,7 章中的 Ising 模型中 的 Gibbs 分布,就是这里的 Gibbs 场的 概率分布的 特殊情形,即 2=L,自旋空间为 }1,1{? 的情形,其 对应的位势是紧邻的 ( 而把 外场看成格点对自己的作用的和 ),Ising 模型中 Gibbs 分布中有一 个倒温度参数,这是为了调节用的,这里无妨 省 去,
再则,由于在不同的邻域中 的科学家有 他们习惯用的记号,所以,在本节中的记号与第
6,7 章中用的有些不同,请读者见谅由 于 记号改动所带来的不便,我们列出其对比,
pba
bpahx
,
,,7,6
本节章第倒温度配称分布组态
t
tX
x
随机场
,
在第 6,7 章中,由于 Ising 模型较为简单,我们只提出了 Gibbs 分布 的概念,而并未提出随机场 的 概念,事实上,Gibbs 分布是 Gibbs 随机场的 不变 分布 （ 有的文献中 Gibbs 分布与
Gibbs 随机场的术语是混用的）,因此 在第 6,7 章中,我们如很多书上那样用 hx,表示组态,
而在本节中,我们强调的是随机场,它用 hx,表示更为合适,所以我们就与很多书上不一样,
改用 ba,表示组态,
定理 9.1 １ (Hammersley - Clifford 基本 定理 )
对于给定的相邻系统?, Markov 随机场与 邻位势 Gibbs 场 是 一样的,】
这个定理的意义在于,用 Markov 随机场来描述图象是一个直观认识,而用 邻位势 Gibbs 场 则 可 归结为其能量函数,后者有时可以用简单的参数来表达 ( 例如相互作用具有平移不变性的时候 ),那时就可用几个参数来描述图象,从而可以达到信息压缩的目的,于
233
是,此定理将直观认识转化为参数表示,以便于数学处理,
１,３ 图象处理的随机过程方法的思路 原则 概述
正如前面所述，一幅黑白图象 是 某个组态空间中的一个组态,在图像传输过程中,需要对于灰度图的信息进行处理与压缩,对于接收到的信息还需要恢复（ Retrieving ）为原图像．常见的方法有,
第 1 种方法 ( 用能量函数的 Gibbs 分布对图像建模 )
按黑白图像的相对灰度 ( 即每个格点的灰度除以图像的总灰度和 )大小,把图像看成格点 集 G 上的组态空间 GLS },,1,0{ L= 中的一个组态，或者看成 G 上的一个概率分布,由于自然景物等的图象往往具有一些局部相联系 ( 空间不变性 ) 的结构,这些局部结构可以抽象地归结为组态间按某种相邻系统的能量函数,
从而可以把一个图像用一个能量函数取到最小 值 的组态，也就是用 Gibbs 分布 达到 最大 值 处的组态，从而可以通过能量函数或 Gibbs 分布来对图像建模,这个 方法 的要害处是：如何 由 图像给出其能量函数,能 量函数表示图像能起一定的信息压缩作用．而反过来由能量函数恢复图像的问题,则 可以用模拟退火或其它
Markov Monte Carlo 方法运作,
第 2 种方法 ( 用低阶 Markov 场的样本对图像建模 )
将黑白图像看成随机场的一个样本,由定理９,１１知道任意随机场都可以看成 Markov 场,而图像的局部结构又说明,用较低阶相邻的 Markov 场的样本来描述图象更为合理,
第 3 种方法 (用迭代随机映射的不动点对图像建模 )
对于一个黑白轮廓图像，构造 与之相系的 有限个压缩映射，使得由此 压缩映射组 给出 的 一个迭 代映射的不动点近似地是此轮廓图．而对于一个黑白灰度图，构造与之相系的一个随机地取有限个压缩映射的随机映射，使此图像近似地是此 随机 映射的不变分布．（ 这 在第３节中将较多地介绍）,这种方法实际上是处理图像的分形方法的另一角度的发展,
一般地，彩色图可以看成红，黄，兰三张灰度图的叠置,
收到了不清晰的图像,也就是 受 到了干扰的图像,就需要对图像进行滤波,称为 图像的清洗,其中常见的一种清洗方法为 Bayes 方法,即把未清洗的图像看成先验分布,而寻找转移机制,使清洗后的恢复图象是此先验分布的 Bayes 分布 (后验分布 ),图像受到污染的典型情形,例如,被 测量的图象象素 xa 受到检测与录制系统的干扰,使 实际 测到的象素为
∑
∈
+=
)(
)(
xy
xxyx whg ab,
其中 g 为一个已知函数,yh 为已知的通道系数,而噪声 xw 的统计规律为已知,这时要恢复图象象素 xa
就需要进行滤波,
此外,图象边 缘信息也可能受到干扰而产生变形,也需要用边缘检测来恢复图象的信息,因此有时应在图象的先验信息中加进,边 -随机场,的先验分布,
有时 在处理图像前可能 还需要对图像进行 有部分相重叠的 分割 (segmentation),
１,４ Gibbs 分布的 样本的 Gibbs 采样法 (Gibbs Sampler)
要对 Gibbs 场作 一些统计 估计 （例如参数估计）,就 先 要得到 Gibbs 分布的配分常数,但是由于组态空间 GLS },,1{ L= 的元素数通常非常 庞大,例如 可以 有 61010,这就使 直接 计算
234
配分常数 Z 时 的求和项太多，从而 常会出现大量被求和的项的值小于计算机误差的精度，
以致使计算 在实际上 失去 可能,用 Markov 链 Monte Carlo 方法 (M C MC) 可以解决这个困难，
即先 近似地得到 Gibbs 分布,再 进行参数 估计,所以 MCMC 是 一个重要的计算手段,
M C MC 方法就是构造以 Gibbs 分布为可逆分布的 Markov 链的转移矩阵,因为这里的随机场的 Gibbs 分布远比 Ising 模型的 Gibbs 分布 复杂,我们只构造离散时间的 Markov 转移矩阵,而且还希望如 Glauber 动力学那样,每次只在一个格点上改变 （ 称为 按 Gibbs 方式 转移 ）,
最后还要保证由这样得到的转移矩阵的各行都收敛于 Gibbs 分布,
为了使构造更为透明,我们先对任给的格点子集 J 定义 如下的,部分转移,
=?
)(0
),()(1),( )( \
其它外坐标相同在 JeZp JSJUH
JJ
baaba ab,(9,4)
其中 JSJ \ab 表示,把组态 a 中在 J 处的坐标换成组态 b 中的相应坐标,而 )(aJZ 是一个归一化常数 (,局部 的,配分常数 ),
∑?=
J
JSJUH
J eZ
b
aba )( \)(,(9,5)
这种,部分转移,表示 组态间只在 位置 J 外的 坐标 相同时才可能 转移,
在此基础上 我们 定义转移矩阵 如下,记 G 中的全部格点的 一个给定的排序为
)|(|),,,{ 21||1 NNGxx G =L,任给组态 S∈ba,,定义
P J SJp ∈= baba,)),((,=P ∏
=
||
1
G
k
P }(
kx
,(9,6)
则 P是一个转移矩阵,注意由于 P }(
kx
中的转移只在一个格点 ix 上进行,由 P构造的 Markov
链 的一次转移,可以视为 || G 次连续转移的结果,其中 每次转移只在一个格点上进行,而且要做遍所有 的格点,
下面的 遍历 收敛 定理是 Ma rkov 链 的 Monte Carlo 方法的 理论 基础,而 Markov 链 Monte
Carlo 方法主要用于获得 Gibbs 分布 的样本,以便进一步得到 Gibbs 场的各种统计平均量,
对于任意格点 Gx ∈,任意组态 S∈a,记
}(min)( }\{}\{ xGxUxxGxU HmH x ahag h
DD
==,(9,7)
那么
),(}{ baxp ∑
∈
=
Gy
mH
mH
xyGyU
xxGxU
e
e
])([
])([
}{\
}{\
ab
ab
=
≥=≥ eGeGGe x
xGxGUU HH
||
1
||
1
||
):)()(max( }{\}{\
d
gbgb
,
其中 xd 为 能量函数 UH 在 格点 x 上的振幅,?为 UH 的 振幅,其 确切 定义为
235
|)()(|max
}{\}{\
bad ba UUx HH
xGxG
= =,xGx d∈?=? max,(9,8)
由此 推出
D≤?≤ epGPC
xx 1),(min||1)( }{,}{ baba,(9,9)
从而 得到
||
}{}{ )1()()()( 1
G
xx ePCPCPC M
≤= L ||1 Ge≤,(9,10)
定理 9,１２ （ 遍历 收敛定理） 以上定义的 转移矩阵 P的 Dobrushin 收缩系数满足
||1)( GePC≤ (?为
UH 的 振幅 ),
记以 P为转移矩阵的 Markov 链为 nx,则它以 Gibbs 分布 p 为可逆分布,而且有
P pTnn 1 →? ∞→,
证明 显见 ),()(),()( abbpbaap JJ pp =,由 P Sp ∈= baba,)),(( 的定义立得
),()(),()( abbpbaap pp =,因此 p 为可逆分布,又由 1)( <PC,用 Dobrushin 定理可知极限成立,
[ 注 ] G - 格点的排序 ),,{ ||1 Gxx L 并不是实质的,我们可以用一个 G 上取值全正的预选概率分布 )(xg,( 即 1)(,0)( => ∑
∈Gx
xgxg 例如均匀分布 ) 的 随机数作 随机选取来代替,这时转移矩阵 P应作相应的改动,代替 ),(}{ ba
kx
p 的是,以 )(xg 的概率对只在 x 坐标不同的组态间进行转移,即令
P (?= ~P ||) G,
~P
Sp ∈= baba,
~
)),((,
其中
==
)(0
)(),()(),( }\{}\{}{~
其它情形若 xGxGxpxgp bababa
,
这时 我们仍有
P pTnn 1 →? ∞→,
１,５ Gibbs 分布的模拟退火
本段 讨论 求 得达到 组态空间 上 能量函数最小值的组态 的概率方法,就是用能量函数的
Gibbs 分布作模拟退火的方法,考虑依赖于温度 nT 的 Gibbs 分布,以它们作为 不变分布构造各步 为 非时齐 Gibbs 转移 矩阵，由此 确定的 一个 非时齐 M arkov 链 }{ nx,我们可以用它作模拟退火,使得当 n 充分大以后，此非时齐的 Markov 链 nx 以 大概率落在 能量函数最小的组态
236
所组成的 集合 中,
更 具体地说,对于 如下的 依赖于温度 nT 的 邻位势 Gibbs 分布
)(1)( 1
)(
a
ap Un
HT
n
n e
Z
=,∑
∈
=
S
HT
n
U
neZ
a
a )(1
,
相应地定义
=?
)(0
),()(1),(
)(1
,
)(
\
其它外坐标相同在 JeZp JSJUn
HT
Jn
n
J
baaba
ab
,(9,11)
∑?=
J
JSJU
n
HT
Jn eZ
b
ab
a
)(1
,
\)(
,
再令
P ()( }{
D
=nx Snxp ∈= baba,)( }{ )),((,P =)(n ∏
∈Gx
P )( }{nx,(9,12)
则 P )(n 是一个转移矩阵,我们 将 以它为 非时齐转移矩阵 列 所确定的非时齐 M arkov 链记为
}{ nx,
定理９,１３ （ Gibbs 分布 的模拟退火的收敛性定理 ）
如果温度 nT 满足 如下的 Geman - Geman 条件,对 00 ↓< nT,存在 0n 使 0nn ≥ 后有
n
GT
n ln
||?≥,(9,13)
其中 D 是能量函数 )(aUH 的振幅,
)(min)(max aaD aaD HH SUS ∈∈?=,
那么在 任意 初 始 分布 0m 下,以 P )(n 为第 n 步转移矩阵列的非时齐的 Markov 链 }{ nx,在时刻 n 的绝对分布 nm 有极限
)(∞→ pm
n,
其中 )(∞p 为集合 min})(:{ =∈ aa UHS 上的离散均匀分 布,从而 还 有
1))(min)((lim == ∈∞→ ax a USnUn HHP,
证明 我们验证 Dobrushin - Isaacson - Madsen 定理的条件成立,条件 (A.1) 的验证与第
8 章 3,2 段的模拟退火定理中的证明完全一样,下面验证 Dobrushin 条件 (A,2) ’,这里定理９,１２中 的不等式变为
237
≤ nT
G
n ePC
||
)( 1)(,
其中?为能量函数 UH 的振幅,于是对于 )()1()( njj PPP L+,我们有
≤+ )( )()1()( njj PPPC L ≤+ )()()( )()1()( njj PCPCPC L )1(
||||
nT
GG
jk
e
=
∏,
由 无穷乘积的理论,当 ∞→n 时,上面右方趋于 0 的充要条件为级数 ∑
n
T
G
ne
||
是发散的,
但是由 Gemam - Gemam 条件,我们有
∑
n
T
G
ne
||
∑?≥
n
ne ln ∑ ∞==
n n
1,
所以 Dobrushin 条件 (A,2)’成立,
［注 1 ］ Gibbs 场的模拟退火 在 图象的清污 ( 图象的滤波 )中的应用
设图象为组态 a,而在传输过程这受到了污染,使得接受方观测到的图象为组态 b,需要用数学 方法进行 滤波,以 便 恢复图象的清晰度,我们把在观测为组态 b 的条件下,状态 a 的分布,记为 )|( baP
(就是后验分布 Posteriori distribution,即 Bayes 分布 ),一般它也可用写成 Gibbs 分布的形式,
)|(
)(
1)|( ba
bba
He
ZP
=,(9,14)
其中 )|( baH 称为 条件能量函数,典型的例子如
)|( baH ∑ ∑
∈
+?=
)(
2
2
2 ))()((
2
1))()((
yx x
xxyx basaaq,
其中第一项强调了图形的光滑性,第二项为保真度,滤 波的目的 就是要找组态
^a
,使条件能量函数达到最小,
)|(min)|( ^ baba a HH =,
这个最佳组态
^a
可以 通过 模拟退火方法求得,
又例如,在边缘检测中,b 就是离散化点上的象素,而 a 就是,边随机链,的一个样值,与滤波问题的不同的是,此时 的 后验图形只是先验 图形的一部分,即边的位置,在复杂的图形处理技术中,例如,X 断层摄影术的处理,情况就远 为复杂,适当的利用分割技术可以得到较好的恢复效果,
[ 注 2 ] 图象的处理有许多数学方法,Gibb s 分布方法只是其中的一种,另一种更为常用的方法是,把图象的一些特征参数作为观察量 ( 这同时也达到了图象信息压缩的目的,用它们来估计 Ma rkov 链 的状态,
就是 重建图象,这就是隐 Markov 方法 (HMM),一般地说,如果特征参数取得合适,HMM 方法是十分有效的,
例如,纹理图象分析方法 (Texture Analysis) 提取特征参数,其它的方法 还 有,下节将介绍的随机函数迭
238
代系统方法,除此以 外,还有小波分析等非概率方法,
[ 注 3 ] 当格点集 G 为无限集时,组态的数目就不再是可数的了,此时就应归入连续状态的 Markov
场的范畴,在统计物理中应用的模型,就是这种模型,有时自旋空间 S 也可以是连续的,例如,圆周 或其它集合,这时就可能存在多个不变分布,称为有相变,因为每个分布代表了统计物理中质点的一个宏观
“相,,统计物理中也正是对这种相变的存在性更感兴趣,
[ 注 4 ] 若每次发展的修改 的 格点集合 J 不是单点集,且各次的修改集可 以 不同,那么,此时的 Markov
链未必有遍历定理,这 时并 不能保证此 Markov 链 nx 的时间平均收敛到不变分布,因此,即使 n 很大,也不能认为 nx 的样本近似为不变 分布的样本,
2 时间离散状态连续的 Markov 链
2,1 概率空间再访
设 }:{ 基本事件?=? ww,F 是样本空间? 的某些子集组成的集合类 (事件体 ),满足
对于任意 ∈nAA,F )1( ≥n,恒有 IU
∞
=
∞
=
∈?=
11
,,
n
n
n
n AAAA W
D
F,
( 这样的 事件体 F 代表了可以通过理论推算,或实际测量能够得到概率的全体 随机 事件 ),在概率理论中 F 称为? 上的一个?s 代数,(?,F ) 称为 概率空间,于是概率 P 就是定义在
s 代数 F 上的一个非负函数,
定义９,１４ （ Borel 集与 Borel 函数）
设 B 是 dR 的某些子集组成的一个?s 代数，且满足以下条件,
(1) dR 中的 任意开集 ∈G B,
(2) 任意一个有 dR 的某些子集组成的?s 代数 G,只要它含有一切开集,则就有
B? G,
那么 B 称为 dR 上 Borel 集合类,其中的集合称为 Borel 集,即 Borel 集合类是以开集为元素的最小?s 代数,( 以前在本书中的所谓,常见的,集合,均指 Borel 集 ),1R 上的 Borel
集称 为一维 Borel 集,
再设 M 为 dR 上的一个满足下述条件的函数类,
(1) 对极限封闭,若 ∈nf M 且 )()( xfxf n →,则有 ∈f M,
(2) 任意连续函数 ∈f M,
(3) 若 dR 上一个函数类 L,只要对极限封闭,而且包含全体连续函数,则有 M? L,
239
那么 M 称为 dR 上 Borel 函数类,其中的函数称为 dR 上的 Borel 函数,( 以前在本书中的所谓,常见的,函数,均指 Borel 函数 ),
再则,随机变量的概念是依赖于 事件体 (?s 代数 ) F 的,严格地应称为 F - 随机变量,
定义９,１５ 样本空间? 上的实值函数 x 称为 F - 可测的,如果对于任意实数 a,恒有 ∈≤ })(:{ awxw F,于是,x 为 F - 随机变量与 x 为 F - 可测的是一样的,而随机变量的分布则依赖于概率 P 的給法,
分布函数为 )( xF 的 随机变量 x 的函数 )(xg 的 数学期望 是如下的 Stieltjes 积分
∫= )()()( xdFxgEg x ( def= ∫ )()( dxFxg )
( 括号中的 )(dxF 就是 )(xdF,前一种记号便于推广到多变量,或条件分布函数 的情形 ),
而
)|)(()|()|( ],( hxhxh xIExPxF?∞
=≤= = hhx =?∞ = yx yIE )]|)(([ ],(
称为 x 对于 h 的 条件分布函数,类似地,我们有
∫= )|()(]|)([ hhx dyFyggE,
( 注意,右方的积分是直观地理解的,如果想要一个清晰的证明,在数学上还必须要作一系列,后台操作,),
2,2 时间离散状态连续的 Markov 链
定义９．１６ 实值随机序列 }0:{ ≥nnx 称 为（时间离散状态连续的） Markov 链,如果对于任意 0≥n,任意实数 10,,,,?nxxxy L,有
====≤+ ),,,|( 00111 xxxyP nnnn xxxx L )|( 1 xyP nn =≤+ xx,(9,15)
(这正说明条件分布函数有 Markov 性）．在 ),,,( 01 xxx Lnn+ 有联合密度时,(9.15)等价于 1+nx
的条件分布密度的如下的等式
====
+
),,,|( 0011
1
xxxyp nnn
n
xxxx L )|(
1
xyp n
n
=
+
xx,(9.15)’
与第５章中类似地，这里的 Markov 性可以有如下的多种等价性描述,
Markov 性的等价性质 1
对于任意实 Borel 集 10,,,?nLLL L,以及随机事件 },,{ 1100 ∈∈= nnA LxLx L,恒有
)(),( 11 xPxAP nnnn =∈==∈ ++ xLxxLx,(9,16)
240
当存在联合密度时,(9.16)式也有其相应 的条件分布密度形式,
Markov 性的等价性质 2
对于随机序列 }{ nx 在大于时刻 n 所确定的随机事件 },,{ 11 knknnnB ++++ ∈∈= LxLx L
及随机事件 },,{ 1100 ∈∈= nnA LxLx L,有
)(),( xBPxABP nn === xx (9,17)
或
)()|()( xAPxBPBAP nn === xx,(9,17)’
（ 直观 证明 我们只在一切条件分布密度都存在情形证明 (9.17)式,先看 2=k 情形,由条件概率的性质,推广了的全概率公式与 Markov 性的定义,我们有
(),( PxABP n ==x ),|,2211 Axnnnnn =∈∈ ++++ xLxLx
（直观地） ),,|( 1122 AxP nnnnn =∈∈= ++++ xLxLx ),|( 11 AxP nnn =∈ ++ xLx
dyyfAxyP nnnnn
n
)(),,|( 1122
1 ++++Λ
==Λ∈= ∫
+
xxx )|( 11 xP nnn =Λ∈ ++ xx
dyyfyP nnnn
n
)()|( 1122
1 ++++Λ
=Λ∈= ∫
+
xx )|( 11 xP nnn =Λ∈ ++ xx,
其中 1+nf 是随机变量 1+nx 的分布密度,由于上式右方与随机事件 A 无关,可见等式左方应该等于
)( xBP n =x,对于一般的 k,只需做归纳法,）,
Markov 性的等 价性质 3
在已知“现在”的条件下，“将来”与“过去” 是条件独立的 (此即 (9,17)’)
Markov 性的等价性质 4
对随机序列 }{ nx 及 0,1 ≥≥? nm 及,任意 Borel 集 L 及实数 10,,,?nxxx L,都有
====∈+ ),,,|( 0011 xxxP nnnmn xxxLx L )|( xP nmn =∈+ xLx,(9,18)
（ 直观 证明 设 一切条件分 布密度都存在,m=1 时即为 Markov 性的定义,对 m 作归纳法,假定 m
时 (9.18)正确,今证 m+1 时它也正确,
),,,|( 00111 xxxP nnnmn ===∈++ xxxLx L
dyyfxxxyP nnnnnmn )(),,,|,( 1001111 ++++ ====Λ∈= ∫ xxxxx L
),,,,|( 001111 xxxyP nnnnmn ====Λ∈=+++∫ xxxxx L
dyxxxyf nnnn ),,,|( 00111 ===+ xxx L
利用归纳法假设及 Markov 性的定义,上式简化为
)|( 11 yP nmn =∈= +++∫ xLx dyxyf nn )|(1 =+ x,
241
注意上式 右方与 },,{ 0011 xxnn == xx L 无关,可知 (9.18)式对 1+m 也成立）,
Markov 性的等价性质５
对于任意有界 Borel 函数 g,即任意实数 10,,,?nxxx L,均有
))((),,,)(( 1100 xgExxgE nnnnn ==== +?+ xxxxxx L１,(9,19)
(在 )()( xIxg Λ= 时,(9,19)就成为 (9,15),然而，从 (9,15)推导 (9,19) 需要测度论的基本知识,本书从略．这里我们粗略地指出其朴素思想,就是?+?= ggg,其中 }0{ >+ = ggIg,而 ±g 们可用有限个示性函数的线性组合来近似 ),
注 等价性质 1,2,3,5 诸款都有相应于等价性质 4 的形式，读者可自己列出它们,
Markov 性的等价性质 6 （最一般的形式）
对于任意实 Borel 集 Λ,及只由随机序列 }{ mx 在时刻 n 及其后的信息所决定的随机变量 nh,恒有
)(),,( 11 xPxxP nnnnnn =Λ∈===Λ∈ xhxxh L (9,20)
或
)()),,,( 0011 xExxxE nnnnnn ===== xhxxxh L,(9,21)
注 这个等价条件的内涵十分丰富,其等价性的证明在测度论中非常典型，本书从略,
与离散状态的 Markov 链的转移矩阵类似,在连续状态 情形,我们有概率转移核,
2,3 概率转移核
定义９．１７ 记
=),(),( Lxp mn )(),|( nmxP nm ≥=∈ xLx,(9,22)
这是一个条件分布 (上式右方就是 )|)(( xIE nm =xxL ),
∈==
)(0
)(1)(),(),(
L
LdL D
L x
xxxp nn,
特别地,),()1,( Lxp nn + 称为时刻 n (到时刻 1+n )的 概率 转移核,而 ),(),( Lxp knn + 就称为从
n 出发的 k 步概率转移核,
特别地,若存在 )(),,(),( mnyxp mn <,使
∫= LL),(),( xp mn dyyxp mn ),(),(,(9,23)
那么 ),(.),( yxp mn 就是,在 xn =x 的条件下,mx 的条件分布密度,称为 条件转移密度,
242
（我们在本书中如无特别声明,则恒假定 条件转移密度 存在．）它满足以下性质,
(P.1) )(1),(0 ),( mnyxp mn <≤≤,∫ 1),(),( =dyyxp mn,
(P.2) ( Chapman - Kolmogorov 方程 ) 对于 nm ≥≥? 有
=),(),( yxp ln ),(),( zxp mn∫ dzyzp lm ),(),(,(9,24)
有时也简写成 乘积核 的形式,即记成
=),( lnp?),( mnp ),( lmp,
( (P.1)是显然的,(P,2)的直观含义十分清楚,它就是条件概率的全概率公式 ),
(P,2)的证明 由 Markov 性及条件期望的性质
)|)|)((( nmlIEE xxxL )|),|)((( nnmlIEE xxxxL= )|)(( nlIE xxL=,
从而
),(),( zxp mn∫∫
Λ
dzdyyzp lm ),(),( )|),(( ),( xdyypE nmlm == ∫
Λ
xx
)|)|)((( xIEE nml == xxxL )|)(( xIE nl == xxL dyyxp ln ),(),(∫
Λ
=,
（注 当不存在条件分布密度时,Chapman - Kolmogorov 方程应写为
=),(),( Lxp ln ),(),( dyxp mn∫ ),(),( Lyp lm ),(9,25)
其中积分是 Stieljies 积分）
2,4 时齐的连续状态 Markov 链
定义９．１８ 连续状态的 Markov 链称 为 时齐的,如果其概率转移核 )()1,( L，xp nn + 与
n 无关 (即 1 步转移概率与出发时刻 n 无关,我们把它改记为 ),( Lxp,
在时齐且有条件转移密度的情形,我们简记
),(),( )1,0( yxpyxp
def
= ),()1,( yxp nn +=,),;,0(),()( ymxpyxp m?= ),(),( yxp mnn +=
这时 Chapman - Kolmogorov 方程就简化为
∫=+ dzyzpzxpyxp nmmn ),(),(),( )()()(,(9,26)
在本书中除非特别声明，所考虑的连续状态 Markov 链总是 时 齐的,
(注 在不存在 条件 转移密度时,由 (9,24) 得到
=+ ),(),( Lxp mnn?+ )1,( nnp（?++ )2,1( nnp ),(),( )(),1( LL xpxp mmnmn =? +?+ ）L
它也不依赖出发时刻 n,它是 Markov 链 }0,{ ≥nnx 的 m 步转移核,其中
∫= ),(),(),()2( dyxpypxp LL,
∫∫ == ),(),(),(),(),( )1()1()( dyxpypdyxpypxp mmm LLL,)
在本 书中除非特别声明，所考虑的连续状态 Markov 链总是 时 齐的,
243
连续状态的 Markov 链的条件转移密度 (或不存在密度时的转移核 ),是刻画此链的基本统计特征量,一个 Markov 链的统计性质可由其初始状态 0x （或其初始分布
)()( 00 xPx == xp ）以及其条件转移密度 (或不存在密度时的转移核 )完全确定,这就是下面的定理
定理 9,１９ （ Markov链的有限维分布密度） 若连续状态的时齐 Markov链 }0,{ ≥nnx
的初始状态为 x=0x,其条件转移密度为 ),( yxp,则
(1) 在条件 x=0x 下 ),,( 1 nxx L 的条件联合分布密度为
),(),(),( 1211 nn xxpxxpxxp?L,
如果 0x 具有分布密度 )( xj,那么,),,( 0 nxx L 的联合分布密度为
)( xj ),(),(),( 1211 nn xxpxxpxxp?L,
(2) 对于有界 n 元 Borel 函数 g 有
== )|),,([ 01 xgE n xxx L
nnnn dxdxdxxxgxxpxxpxxpn LLLL 2111211 ),,(),(),(),(12?ΛΛΛ∫ ∫ ∫,
（证明提示 (1),对于 n 用归纳法,(2) 对于 )()( 11 nn xgxgg L= 可用归纳法证明,而对于一般的 g,则需要用测度论中的典 型逼近方法,这里从略,(其粗略思想为,用示性函数的线性组合近似
}0{ >
+ =
ggIg,,再用形如 )()( 11 nxIxI nΛΛ L 的函数的线形组合,去近似一般类型的示性函数
),,( 1 nxxI LΛ,其中 L 是任意 n 维 Borel 集 ),
注 在转移核不存在密度且 0x 没有分布密度 的情形,),,( 0 nxx L 的联合分布为
∫∫
ΛΛ
=Λ∈Λ∈
n
nnP LL
0
),,( 00 xx )(0 dxF ),(),(),( 1211 nn dxxpdxxpdxxp?L,
其中 )(0 xF 是 0x 的分布函数,),( Λxp 是 Markov 链的转移核,而
== )|),,([ 01 xgE n xxx L
),,(),(),(),( 11211
12 nnn
xxgdxxpdxxpdxxp
n
LLL?ΛΛΛ∫ ∫ ∫,
定理９,２０ （连续状态的 Markov 链的绝对概率 ）
设 )|( xyf n 是在条件 x=0x 下,随机变量 nx 的条件分布密度．那么
244
dzzxpzyfxyf nmmn ),()|()|( )( ）（∫=+,(9,27)
证明 应用全概率公式,
2,5 例
例９．２１ ( 非随机的迭代映射 ) 设随机 变量 序列 nx 满足迭代关系
xf nn ==+ 01 ),( xxx,
这时 nx 没有分布密度,事实上它是 一个常数,即由 x 只能转移到 )(xf,于是,对于任意实
Borel 集 L 可以直接计算得到
),,,|( 00111 xxxP nnnn ===∈+ xxxLx L
),,,|)(( 0011 xxxxfP nnn ===∈= xxxL L })({ L∈= xfI,
也就是 其转移核为
)(),( )(1 xIxp f LL?=,
其中 )(1 L?f 是 L 在 f 下的原象集,即 })(:{)(1 LL ∈=? yfyf,
例９．２２ ( 带随机噪声的迭代映射所生成的连续状态的 Markov 链 )
设 }{ nw 为独立序列，且 nw 的分布 密度 为 )( yf n,),( yxg 是 2R 上 的 Borel 函数．又设随机变量序列 nx 满足迭代关系
),(1 nnn wg xx =+,
其中 0x 与 }{ nw 独立,此时 nw,),( nwxg 都与 },,,{ 10 nxxx L 独立．故而
),,,|( 00111 xxxyP nnnn ===≤+ xxxx L
),,,|),(( 0011 xxxywxgP nnnn ===≤= xxx L
∫ ≤=≤= dzzfywxgP nyzxgn )()),(( ),(,
上式右方与 01,,xxn L? 无关,从而等于 )|( 1 xyP nn =≤+ xx,可见 nx 是连续状态的 Markov
链,而且,当且只当 )( yf n 不依赖 n 时,这个连续状态的 Markov 链才是时齐的,
2,6 寻找 dR 上可微函数 )(xf 的最小值的位置的模拟退火算法
先考虑 1=d 情形,设 nw 为独立同分布的 )1,0(N 随机变量序列，又设 函数 0)( ≥xf
245
且二阶光滑，再设 )(xf 取最小值的点全体组成集合 0S,我们考虑如下的随机迭代模型
001 ),0()()(' xhfwhf nnnnnn =>++=+ VVeVVV,
它是例９．２ ２的特殊情形，而且在 xn =V 的条件下,转移到 yn =+1V 的 条件转移 密度为
=+ ),()1,( yxp nn
2]
)(
)('[
2
1
)(2
1 xf xhfxy
n
ne
xf
e
pe
,
即它是正态分布 ))(),('( 22 xfxhfxN ne+,其中 x 作为参数,于是 nV 是时间离散的非时齐的状态连续的 Markov 链,这里的非时齐性是由于“小参数” ne 依赖于时刻 n 而引起的,
对于 1>d 情形,则相应的随机迭代模型为
001 ),0()()( xhfwfh nnnnnn =>+?+=+ VVeVVV,(9,28)
其中 nV,nw 都是随机向量,而 f 为多元函数 ( 注意如果 0Sn ∈V,则有 0)( =? nf V,此时随机迭代模型就变成
001 ),0()( xhfw nnnnn =>+=+ VVeVV,
即只有随机小扰动项 )( nnn fw Ve 起推动 nV 运动的作用,由此可知一旦当此 Markov 链的轨道 nV 到达了 f 的某个稳定点以后,再到它的给定的某个小邻域以外的可能性就较小，而且当此稳定点是 f 的最小值点时，这种可 能性就相对地更小,我们再注意,当 此 Markov 链的轨道 nV 在未到达 f 的某个稳定点之前，在（９．２８）中的项 )( nfh V? 起的主要作用将使
nV 往 f 的稳定点发展．综合此两个朴素思想的启示,我们可以用此 Markov 链的发展来“近似”地找到 f 的最小值的位置，这正是连续 状态空间变量的模拟退火思想的要点,当然用这种方法还是有失败的可能的．然而，这样做最终失败的可能性较小,如果我们能忍受小概率失败的失误而放弃,而只求以较大的概率能达到整体最小值附近,那么这是一种可行的妥协方法 ( 这种加入随机小扰动项 )( nnn fw Ve 的方法，能避免用确定性方法 ( 例如梯度方法 )
求最小值位置时,陷入局部极小处的失误 ),把 ne 看成温度后,（９．２８）提供了连续状态空间的一种模拟退火方法,即利用增加小扰动项的随机迭代模 型作模拟退火,如果适当地选取退火的温度 ne 趋于 0 的速度（例如,如 Geman - Gemam 那样选取 1)(ln?= nn ge,而 g 又足够大）,那么,在理论上可以证明,对于任意 0>d 有
1))(( 0 →?∈ ∞→nn SBP dV,
246
其中 )( 0SBd 指 0S 中所有点的δ邻域的并 ( 这 种极限关系称为 nV 按 概率收敛到集合 0S,与有限状态空间的情形看起来不同,在 有限状态空间 情形 是收敛到最小值集合,而现在时收敛到最小值集合的一个邻域,如果我们认为有限集中的点都 是孤立的,那么这两者是一致的 ),
也就是说,对随机差分方程作迭代,经过足够多次以后,就以近似于 1 的概率,保证 nV 会到达 f 取最小值的集合的较近的范围,这就 是用随机迭代模型作模拟退火的基本思想，这里温度 ne 依赖 n 的阶 1)(ln?= nn ge,对于概率收敛性的成立是至关重要的,
注意求函数 )(xg 零点的问题,可以化为求 2)(xg 的最小值问题,因而也可以应用随机迭代模型的模拟退火方法,因此这个数学模型在广泛的领域中有应用前景,
[ 注 ] 模拟退火方法中的人工随机干扰 nw,是 独 立同分布随机 序列,称为 白噪声,模拟退火的进程 一般较 慢,对此 人们采取了 多种加速收敛措施,其中常用的有 TINA ( Time
Invariant Noise Annealin g,用时间不变的噪 声 作 模拟退火 ) 算法,它的基本思想是,采用,有色噪声,作为逃逸出局部极值 的,驱动力,,即在目标函数值大的时候,使干扰大些,
而当函数值接近最小值时,就把干扰取得小些,特别是在已知最小值的情形 ( 例如求解非线性联立方程组时,只要把目标函数取为各方程 的 平方和,那么目标函数的最小值就是 0,
而 达到最小值的变量就是该方程组的解 ),例如,在已知最小值的条件下,对 )(xf 求整体最小值 的位置,也 可以考虑如下随机差分方程,
nnnnn wfff ))(()( *1?+?+=+ xxxx,(9,29)
其中 nw 是独立同分布的随机变量列,*f 是目标函数 )(xf 的最小值,于是 (９,２９）给出了用有色噪声 nn wff ))(( *?x 干扰的 TINA 算法,
我们还可以将 TINA 算法修改为适合于未知最小值的全局优化问题,例如可以把干扰噪声取为 nn wncfg ])ln()))(([ 1?+?+? gx,其中 g 是一个有界增函数,c 是 一个常数,
它取得越接近最小值 *f 越好,而在实际计算中,人们可以根据对 f 的 函数值的认识,不断地调小参数 c,以求加速收敛,
2,7 Dobrushin 不等式,指数遍历性与收敛性
在第 5 章中，我们曾对离散状态的 Markov 链,讨论了当 n →∞时它的极限性 态,在本节中，我们将利用 Dobrushin 不等式，对于连续状态的 Markov 链得到相应的结果．这类结果要求假定转移密度存在，并满足较强的条件，这时不仅可以得到 n 步转移密度当 n →∞ 时的极 限趋于唯一的不变概率分布密度，而且保证误差具有指数型的下降速度 ( 称 之 为指数遍历性 ),这一结论比之于一般 Markov 链的遍历定理要强得多,
命题 9,２３ ( Dobrushin 不等式 )
对于 Borel 函数 f,定义 dxxff |)(||||| ∫=,假定时间离散的连续状态的 Markov 链
247
存在转移密度 ),( yxp,那么对于任意分布密度 )(),( 21 xx rr 有
||||)(||)(),()(),(|| 2121 rrrr?≤ ∫∫ PCdxxxpdxxxp,(9,30)
其中
dzzypzxpPC yx |),(),(|sup21)(,?= ∫ （９．３１）
称为转移密度 ),( yxp 的 Dobrushin 数,】
例９,２４ 设 )(xf 是一个有界连续函数．那么对于满足 )()( yfxf ≠ 的点 ),( yx,
2
)()(),( yfxfyxh +=L 是方程
2
))(( 2xfz
e
2
))(( 2yfz
e
=
的唯一解,且 ),( yxh 也是有界连续函数,当 )()( yfxf < 时我们有
∫ ≤? dzzypzxp |),(),(| ),( yxh∞?∫ p21 dze
xfz
2
))(( 2 ∞
∫+ ),( yxh p21 dze
yfz
2
))(( 2
)(),( xfyxh?
∞?∫
＝ p21 due
u
2
2
∞
∫
+
)(),( yfyxh p2
1 due u22?
2m
∞?∫
≤ p21 due
u
2
2
∞∫+
1m p2
1 due u22? ))(1()(
12 mm Φ?+Φ=,
其中
}|:)(),({|inf,1 yxxfyxhm yx ≠?=
D
,}|:)(),({|sup,2 yxxfyxhm yx ≠?=
D
,
∫
∞?
=
a
a)(F p21 due
u
2
2
,
显见 ∞<<<∞? 21 mm,于是
yxPC,sup2
1)( =
p2
1 2 ))(( 2| xfze∫ dze yfz |2 ))(( 2
2
1≤ 1))](1()([
12 <Φ?+Φ mm,
例９．２５ 设连续状态的 Markov 链的转移密度 ),( yxp 满足条件
0)(),( ≥≥ ygyxp,
那么此转移密度的 Dobrushin 数满足
248
1)(1)( ≤?≤ ∫ dyygPC,
证明 仿照第５章离散状态的 Markov 链的 Dobrushin 不等式的证明可知
=?= ∫ dzzypzxpPC yx |),(),(|sup21)(,
= dzzypzxpzypzxpyx )),(),((sup ),(),(,?∫ ≥ dzzgzxpzypzxpyx ))(),((sup ),(),(,?≤ ∫ ≥
≤ Sup dzzgdzzgzxpyx )(1))(),((,∫∫?=?,
定义９．２６ 具有转移密度 ),( yxp 的连续状态的时齐 Markov 链 nx 称为 指数 L 1遍历 的 ( 或 强指数遍历 的 ),如果存在正整数 0n,分布密度 )( y）（ ∞r 及常数 0,10 ><< cr,只要 0nn ≥,就对任意 x 满足
nn crdyyyxp <? ∞∫ |)(),(| )()( r,
其中 ),()( yxp n 是此 Markov 链的 n 步转移密度,
下面的概念只需要转移核，而不需要假定存在转移密度,
定义９．２７ 连续状态的时齐马氏链 nx 称为 指数遍历 的，如果存在正整数 0n,分布函数 )()( yF ∞ 及常数 0,10 ><< cr,只要 0nn ≥,就对任意 yx,满足
n
n cryFyxP <?≤
∞ |)())((| )(x,
这里 改 用 )(xnx 记初值为 x=0x 的 Markov 链 nx,分布函数 )()( yF ∞ 不依赖于初值 x,称为
nx 的 ( 遍历 ) 极限分布,也称为 nx 的 平稳分布,以
)(∞F 为初分布，就能得到一个连续状态的平稳过程,显见,指数 1L 遍历一定指数遍历,
定理９．２８ （指数 L 1遍历定理）
具有转移密度的连续状态时齐 Markov 链 nx,如果存在 10 ≥n,使其 0n 步转移密度
),()( 0 yxp n 的 Dobrushin 数 1)( )( 0 <nPC,那么此 Markov 链是指数 L 1遍历的,
( 证明 与第 5 章中离散状态情形 类似,从略 ),
例９．２９ 例９．２５中 的 Markov 链,只要 存在 10 ≥n,使其 0n 步转移密度满足
),()( 0 yxp n )( yg≥,且 ∫ > 0)( dyyg,
则它 指数 L1遍历
249
例９．３０ 设 f 是有界函数，而 nw 是独立同分布随机变量序列，且具有严 格正的连续的分布密度 )( xp,那么由如下定义的迭代列 nx 所决定的连续状态 Markov 链是指数 L 1
遍历的
nnn wf +=+ )(1 xx,(9,32)
证明 不妨设 bxfa ≤≤ )(,此 Markov 链的转移密度为
min))((),( ≥?= xfypyxp bza ≤≤ )( zyp? ＞ 0,
由例９,２５和定理９．２８可知此 Markov 链是指数 L1遍历的,
定理９,３１ （指数遍历定理）
设 nw 是独立同分布随机变量序列,),( vuf 是连续函数，且关于 u 单调不减,又
)),(()(,10 nnn wxfxx xxx == +,(9,33)
如果存在 z 0,0n 及 0>d 使
))((sup 0
0
zxP nx ≤x 及 ))((inf 0
0
zxP nx ≥x 都 d≥,
则连续状态的 Markov 链 )(xnx 是指数遍历的,
( 证略 ),
3 随机的迭代函数系统 ( 随机 IFS,随机 Iterative Function System)
3,1 局部相似性的基本想 法
如第２节所示，图象的局部相似性是它的基本特点,最简单的图象是只有黑白两个灰度的轮廓图,由于这种图的轮廓线在局部总有一点相似性 ( 所谓自相似性 ),人们就希望用函数迭代映射的,不动点,,来近似描绘轮廓线图,并进一步把一般的黑白灰度图和彩色图用随机迭代函数映射的不变分布来近似,这样就可以达到在传输过程中压缩图象信息,而在接收后再重建合成人工图象的目的,( 显然这是非常朴素的想法,例如,压缩图象信息的压缩比如果能达到 20,则效果就很不错了,但是在实践中,还必须考虑到与原有的设备能相容这个方 法,即改进压缩比的方法只希望在原有设备上附加一个系统,而不希望完全改建设备,因此再理想的压缩方法也必须与实际情况相协调 ),
图像的纹理分析 (Texture Analysis)是一种重要的图像处理方法,在目标识别等方面有重要应用。图象的纹理是一种结构,它包括：自然纹理表面的灰度，颜色，条纹,以及它们的模拟,要模拟或分析图像纹理，选择合适的模型是关键的一步,纹理模型必须能表述图像元素的空间分布特征,均匀的纹理的最基本的分布特性是它的空间相似性，就是在图像上的不同地方的纹理有着几何上相似的象素组 合,在数学上，可以用小波分析,自相似随机过程,
随机迭代函数系统,Markov 随机场等来描述纹理,由于 Markov 随机场和 Gibbs 随机场是等价的，因而就可以用 Gibbs 随机场模型 (用另一种看法视为随机图,而把图像作为随机图的
250
一个样本 ) 描述纹理．进而可以使用各种方法估计 Gibbs 分布的参数，在此基础上就可以对纹理进行分析和处理,下面介绍的轮廓图的方法,是一种与 Gibbs 分布不同的,近似描绘 无灰度的黑白图的方法，即迭代函数系统方法,
3,2 轮廓图全体组成的距离空间
定义９,３２ （轮廓图的图象窗口中的距离）
只有轮廓而无灰度的黑白图,即轮廓图,其黑色部分就可看成为平面上某个窗口的一个闭集,我们记此窗口为 X,它是 2R 中的有界闭集,例如,]1,0[]1,0[ ×=X,X 中的两个点 ),(),,( 2121 yyyxxx == 之间可定义距离 ),0[,∞→× XXd,例如普通的 Euclid 距离 222211 )()(),( yxyxyxd?+?= ( 也可以是其它的距离,只要满足如下的 3 条距离公理,
(D.1) 0),( ≥yxd,而且等号成立的充要条件为 yx = ;
(D.2) ( 对称性 ) ),(),( xydyxd = ;
(D.3) ( 三角形不等式 ) ),(),(),( zydyxdzxd +≤,
例如,|)||,max(|),( 2211 yxyxyxd=,或者 ||||),( 2211 yxyxyxd?+?= ),
轮廓图的图象空间与图像之间的距离
如上所述,一幅轮廓图就是窗口的一个闭子集,从数学上描述两幅轮廓图的近似程度的最朴素的想法是,在两幅图之间引入距离 ( 也可以引入准距离,相对熵,,这将在第 10 章的隐 Markov 模型与第 15 章的算法中叙述 ),为此,我们先介绍一些关于距离空间的基本知识,由前面的论述,我们知道,窗口 X 可配以距离 d,我们把 ),( dX 作为蓝本,引入距离空间的概念如下,
定义９．３３ 任意一个集合 Y,如果定义了满足 3 条距离公理的映射 ),0[,∞→× YYd,
则我们称 ),( dY 为 距离空间 (metric space),距离空间中有收敛概念,对于 Yxn ∈,如果有 0),( →xxd n,则称 nx 收敛到 Yx ∈,记为 xxn →,
距离空间中的序列 )(},{ Yxx nn ∈ 称为 Cauchy 列,如果对于任给 0>e,存在 N,只要 Nmn >,,就有 e<),( mn xxd ( 显见这个概念是欧氏空间中的 Cauchy 序列的推广,需要注意的是,在距离空间中并没有加法运算,除非它还具有线性空间结构 ),
距离空间 ),( dY 称为 完备的,如果它的任意 Cauchy 列 }{ nx 都有极限,即存在 Yx ∈,
使 0),( →xxd n,( 例如,dR,dR 中带边的超立方体,dR 中带边的有界图形,都是完
251
备的距离空 间 ),
对于距离空间 ),( dY,Y 的子集 F 称为 闭集,如果 F 中的任意收敛序列的极限仍在 F
之中,闭集的余集称为 开集,Y 的子集 K 称为 紧集,如果 K 中的任意序列都有收敛 子序列,
而且此收敛子列的极限仍在 K 之中,
显见紧集必是闭集,在欧氏空间中,紧集与有界闭集是一样的,在上述的窗口空间中,
由于窗口本身是有界的,因而闭集与紧集是一样的,但是,在一般距离空间闭集就未必是紧集,
定义９,３４ （作为距离空间的轮廓图的图象空间）
设 ),( dX 为窗口距离空间,在此窗口空间上的一个轮廓图,就是窗口空间的一个闭集,
记
}:{)( XFFXH?=? 闭集,
那么,)( XH 为所有的轮廓图的图象组成的集合,称为 图象空间,对于任意两个轮廓图
)(,XHBA ∈,可以定义如下三种不同的距离 21,,rrr,
),(sup),(sup),( AydBxdBA ByAx ∈∈
+=r,
其中
),(inf),( yxdBxd By∈= ;
},:{inf),(1 ddd dr ABBABA=?,
其中
}),(:{ dd <∈?= yxdAxyA 使 ;
)],([sup)],([sup),(2 BydBxdBA ByAx ∈∈
∨=r,(9,34)
其中
),max( baba =∨,
不难验证 21,,rrr 都是 )( XH 上的距离,而且它们都彼此等价,即收敛性是一样的 ( 虽然距离的大小并不一样 ),
命题９．２５ 距离空间 )),(( rXH 是完备的,
证明 设 nA 为 )),(( rXH 的 Cauchy 列,则易见对于任意 nn Ax ∈,nx 是 ),( dX 的 Cauchy 列,
故有极限 x,令
})(:{ xAxxA nn →∈?=?,
用与普通数学分析类似的方法,可以证明 A是闭集,即 )( XHA ∈,而且 0),( →AAnr,这就证明了
252
)),(( rXH 的完备性,
完备距离空间的压缩映射与不动点
定义９．２６ 任意距离空间 ),( dX 自身的映射 W,如果存在 1<c,使得对于任意
Xyx ∈,,恒有
),())(),(( yxcdyWxWd ≤,(9,35)
那么,W 称为 压缩映射,c 称为 压缩系数,
例９．２７ (仿射压缩映射 ) 设 W 为如下映射
+
→?
=
2
1
2
1
2221
1211
2
1
a
a
x
x
aa
aa
x
xx W,
假定系数
2221
1211
aa
aa,
2
1
a
a 满足：
1
1,
0
1,
1
0,
0
0 的 W 像都在 X 中,且
=c
2221
1211
aa
aa T
aa
aa
2221
1211 的最大特征值 1<,那么,可以证明 这样的仿射映射是压缩映射,
定理９．２８ ( 不动点定理 ) 对于任意完备距离空间 ),( dX 上的压缩映射 W 有
(1) 存在唯一的不动点 aaWa =)(,,
(2) (迭代近似 ) 记,2 WWW oo = )),(()(2 xWWxW =o,1 nn WWW oo o=+
则对于任意 Xx ∈,恒有 axW n →)(o,
(3) (近似不动点 ) 如果 y 满足 e<)),(( yyWd,则 cayd?≤ 1),( e,
证明 先证明不动点的存在性,令 =nx )(xW no,于是 ),(),( 11 xxdcxxd nnn ≤+,从而
∑
=
+++ →?≤≤
p
k
n
knknnpn xxdc
cxxdxxd
1
11 0),(1),(),(,
所以 }{ nx 是 Cauchy 列,由完备性,它有极限,记为 a,再则
0),()),(())(),(()),(( →++≤ axdxxWdxWaWdaaWd nnnn,
故 aaW =)(,再证明不动点只有一个,假 设 另有不动点 b,则 有
),())(),((),( bacdbWaWdbad ≤=,
253
由 1<c 立得 0),( =bad,即 ba =,这就证明了 (1)和 (2),最后
)),(())(),(())(,(),( 1 ayWdyWyWdyWydayd nnn oooL +++≤?
)),((1 ayWdcc nn oL ++++≤? eee c?≤ 1 e,】
从窗口空间 ),( dX 的压缩映射导出图象空间 )),(( 1rXH 的压缩映射
定理２．２９ ( Collage 拼贴定理 )
若窗口空间 ),( dX 上有有限个压缩系数为 ic 的压缩映射 iW )1( Ni ≤≤,定义
)),(( 1rXH 上的映射
^W
,其象集为所有的 )( NiWi ≤ 的象集之并,
)( XHA ∈?,U
N
i
i AWAW
1
^
)()(
=
=,
那么,
^W
是 )),(( rXH 上的压缩映射,
^W
称为由 },,{ 1 NWW L 导出的 迭代函数系统,它的唯一不动点 0A 称为 迭代 函数 系统
^W
的 吸引子,而且对于任意 )( XHF ∈ 有误差估计 ( 称为拼贴定理 ) ( 设 ic 们中的最大数为 c )
))(,(1 1),( ^101 FWFcAF rr?≤,(9,3 6 )
证明 由 )( XH 上的等价距离 1r 的定义得到
))(),(( ^^1 BWAWr U U
N
i
N
i
ii BWAW
1 1
1 ))(),((
= =
= r
))(),((max 1 BWAW iii r≤ ),(1 BAcr≤,】
迭代函数系统的吸引子的＂自相似＂特性
由 },,{ 1 NWW L 导出的迭代函数系统
^W
的吸引子 0A 满足,
UN
j
ji AAWAWAW
1
00
^
00 )()()(
=
==?,
这说明 吸引子 0A 的 iW 象是 0A 的一部分,再用 kW 们继续地迭置作用下去,则它相当于对
0A 不断地作细分,在国内外工程文献中,称以上性质为,迭代函数系统
^W
的吸引子 0A 关于 },,{ 1 NWW L 有自相似性,事实上，当 kW 们都是仿射压缩映射时，这种说法较为合理,
而在一般情形是一种近似的描述,在实际运作中，迭代函数系统中的映 射常取成仿射压缩
254
映射，这就充分反映了轮廓图的自相似性,
用迭代函数系统压缩图象信息的思路为：通过选取 N 及 },,{ 1 NWW L,用迭代函数系统
^W
的吸引子 ( 在距离 1r 下 ) 去近似标的图象．由此可以达到压缩信息的目的，因为这时轮廓图的信息完全可由 },,{ 1 NWW L 描述,下面的定理说明上述要求不但能够做到,而且可以把 iW 们选取为仿射映射,
定理９,３０ ( Hutchinson 定理 )
任给的轮廓图 )( XHF ∈ 及容许误差 0>e,存在 1>m,及压缩的仿射映射组
},,{ 1 mWW L,使
^W
的吸引子 0A 是 F 的容许 近似,即
er ≤),( 01 AF,
证明 因为 F 是紧集,所以可以用有限个半径为 8e 的圆 mBB,,1 L 覆盖 F,U
m
k
kBF
1=
,我们不妨假定这些圆的中心都在 F 的 8e 邻域中,于是 4
1
e
FBF
m
k
k
=
U,我们可以找压缩系数 21<ic
的仿射映射 iW,使 ii BFW?)(,又由于 iB 的直径小于 4e,故
2)()(
e
FWBFW iii,
从而
4
^ )( eFFW?
,2
^ )( eFWF?
,
因此，由距离 1r 的定义得到
2))(,(
^
1
er ≤FWF
再用 Collage 定理推出
er ≤),( 01 AF,
[注 ] Hutchimson 定理出现于 1981 年,它的主要含义是,可以把轮廓图近似地作为迭代函数系统的吸引子,吸引子具有,自相似性,,也就是工程界广泛认为的,分形性质,,
轮廓图是窗口 X 的一个闭子集，可以用迭代函数系统的吸引子来建模．而窗口上的灰度图的相对灰度可以看成口上的一个概率分布,为此考虑用随机迭代函数系统生成的
Markov 链的不变分布来建模．这就是下一段的理论,
在窗口 ]1,0[]1,0[ ×=X,相对灰度图就相应地考虑为 概率 空间,( X B ))( X 上的一个概率分布,其中 B )( X 为 X 的 Borel 子集全体组成的类,
255
3,3 灰度图与随机迭代函数系统
本段考虑最简单的随机映射 及由它 生成的 Markov 链，并用以建模．这种随机映射就是,对于 有限个映射，以不同的概率从中取一个．其确切的定义为
定义９．３１ 窗口 ]1,0[]1,0[ ×=X 自身的 随机压缩映射 W的定义为,存在 X 自身的压缩映射 )()1(,,mWW L,使
(P W )1(,) 1)( =++== mii pppW L,
直观地,我们可以 将 它 记为 如下的分布 列的形式 （与离散分布的不同处是其取值为映射，而不是数）
W
)()1(
1
,,
,,~
m
m
WW
pp
L
L,(9,37)
我们 称
)()1(
1
,,
,,
m
m
WW
pp
L
L 为 随机映射 W的分布,
命题９,３２ （独立同分布的随机映射生成的 Markov 链）
设 LL,,,1 nWW 为独立同分布的随机映射序列,设 ０x 是随机变量．定义
xxx,),( 011 LW= =+1n W )(1 nn x+,(9,38)
则 }{ nx 是时齐的 Markov 链,称为 随机迭代函数系统,其概率转移核 (不存在条件转移密度 )
为
),( Lxp ∑
=
=
m
i
Wi xIp i
1
)( )(1)( L,(9,39)
证明
),,,|( 00111 xxxP nnnn ===∈+ xxxLx L
(P= W ),,,|)( 00111 xxxx nnnn ===∈+ xxxL L
(P= W ))(1 L∈+ xn ∑
=
∈=
m
i
i
i xWPp
1
)( ))(( L ∑
=
=
m
i
Wi xIp i
1
)( )(1)( L,
推论９．３３
∑
=
+ ==
m
i
i
inn xWfpxfE
1
)(
1 ))((]|)([ xx,】
下面介绍独立同分布的随机压缩映射生成的 Markov 链不变分布的存在与唯一性,
定理９．３４ ( 遍历性定理 )
随机压缩映射生成的 Markov链 nx 存在唯一的不变分布,设它是某个随机变量 V 的分
256
布函数 )xF（V,那么对于此 Markov 链的转移核 ),( Lxp,分 布函数 VF 有 以下的不变性质,
=)(xFV )(]),(,( ydFxyp V∫?∞,
其中 ],(],(],(),,( 2121 xxxxxx?∞×?∞=?∞=?,VF 称为 不变分布函数,而且,对于任意有界 Borel 函数 f,不管什么初值 0x,恒有遍历定理，即有如下的按概率收敛式
)()()()()(1
1
∞→=?→?∑ ∫
=
nxdFxfEffn
n
i
p
k VVx,
特别地,如果把 nkk ≤≤1}{x 中落入 ],( x?∞ 的频数记为 ],( nn?∞,则有以下的按概率收敛式,
)()(],( ∞→?→∞ nxFnn px V,
这就得到了不变分布 函数 的相合估计如下 (如果将一幅灰度图用一个随机迭代函数系统 建模,
那么,利用 此相合估计 可以在实际计算中,通过随机迭代函数系统重建灰度图 )
=)(
^
xFV nn x],(?∞,(9,40)
证明 令 x== 00 xh,定义
hhhhh,),(),( 12011 Lnn WW == + =+1n W )(1 nh,
那么 }{ nh 与 }{ nx 同分布,由于 }{ nh 是某个概率空间,(? F,)P 上的随机变量序列,对固定某个基本事件 (样本 )?∈w,存在 mii n ≤≤,,1 L１ 使 )()( )()( 1 xWW niin oLo=wh (其中运算 o 表示映射的复合 ),于是对于 此固定的基本事件 w,压缩映射 )()1(,,mWW L 中最大的压缩系数 )1(<c 有
))(),(( whwh npnd + )(( )()( 1 xWWd nii oLo= ))((,)()( 1 xWW pnii +oLo
0))(,( )()( 1 →≤ ++ xWWxdc pnn iin oLo,
这说明 ))({ whn 是 X 中的 Cauchy 列,从而它收敛于某个 (依赖于初值 x 的 ) )(wV x,即存在随机变量
xV (理论上可以严格证明 ∈≤ }{ axV F ),使
1)( = →? ∞→ xnnP Vh,
再则,由于 )(whn 也依赖于初值 x,我们 将 它改记为 )(wh xn,又 因为它是压缩映射的复合,所以仍然是压缩映射,故 有
257
0),())(),(( →≤ yxdcd nynxn whwh (关于 w 是 一致 的 ),
由此得到
))(),(( wVwV yxd ))(),(( whwV xnxd≤ ++ ))(),(( whwh ynxnd ))(),(( wVwh yynd 0→,
这说明极限随机变量 xV 与初值 x 无关,所以我们可以 将 它改记为 V,显见,nh 的分布函数
n
Fh 按分布 收敛 到 V 的分布函数 VF ( 即在 )( xFV 的连续点 x 上有 )()( xFxF
n Vh
→ ),
由 nx 与 nh 同分布得到,nx 的分布函数
n
Fx 也 按分布 收敛到 V 的分布函数 VF,可以证明,这个分布 VF 就是转移核 ),( Lxp 的唯一的不变分布 函数,
(注 此 事实可以证明如下,对于任意有界连续函数 f,有 (f x ()1 fn =+ W )(1 nn x+ )
(Ef W ）＝)(1 nn x+ ∑
=
m
i
n
i
i WEfp
1
)( ))(( x ∑
=
=
m
i
n
i
i WfpE
1
)( ))](([ x
∫∑
=
= )())((
1
)( xdFxWfp
n
m
i
i
i x ∫∑?
=
= ))(()( 1)(
1
xWdFxfp i
m
i
i nx,
但是按定义有 (f x ()1 fn =+ W )(1 nn x+ ),对此 两边取期望便得到
)()(
1
xdFxf
n +∫ x ∫ ∑
=
= ))(()( 1)(
1
xWFpdxf i
m
i
i nx,
令 ∞→n,利用 按分布 收敛的性质 推出
)()( xdFxf V∫ ∫ ∑?
=
= ))(()( 1)(
1
xWFpdxf i
m
i
i V,
由 f 的任意性,再利用（ 9.39） 式便 得到
)( xFV ∑
=
=
m
i
i
i xWFp
1
1)( ))((
V )()(
1
],()( 1)( ydFyIp
m
i
xWi i V∑ ∫
=
∞?= )(]),(,( ydFxyp V∫?∞=,
(把它 写成测度的形式就是
=)(LVF )(),( ydFyp VL∫ ),
这就证明了极限随机变量 V 的分布 VF 是转移核 ),( Lxp 的不变分布 函数,再则,如果此转移核另有不变分布 函数 F,将 以它为分布 函数 的随机变量 记为 J,取 Jx =0,那么,由 迭代定义的随机变量序列
}{ nx 就是一个平稳 序列,与 上面的证明 类似地 仍有 1)( =→ VxnP,于是作为极限分布的 F 应该等于
258
VF,此外,还是由于 1)( =→ VxnP 及 平稳列的性质,还 可以 进一步 证明有不依赖初值的遍历定理 ),】
(对于具有一些抽象积分知识的读者,我们 还 可以用距离空间的方法来证明不变分布 函数 的 存在 性如下,
令 M 为窗口空间 ),( dX 上全体概率测度的集合,即 M 中的 每 一个元素 m 都是 X 上的 Borel 集合类
B )( X 上定义的一个概率,在 M 上可以定义如下的 Hutchinson 距离 r,
∫ ∫?= ∈? |)()()()(|sup),( xdxfxdxfLf nmnmr,
其中 )},(),,(|)()(:|{ XyxyxdyfxffL ∈?≤?=?,可以证明 ( M,r ) 是完备距离空间,对于压缩系数为 ic 的压缩映射 )(iW 们组成的随机映射
)()1(
1
,,
,,
m
m
WW
pp
L
L
,我们定义它在 M 上导出的 一个 映射为 (记为 1?W )
∑
=
=→
m
i
i
i WpW
1
1)(1 )( mmm,
往证 1?W 是压缩的,
= ),( 11 nmr WW ∫ ∫∈? |))(()())(()(|sup 11 xWdxfxWdxfLf nm
∫ ∫∑?=
=
∈ |)())(()())((|sup
)()(
1
xdxWfxdxWfp ii
m
i
iLf nm
∫ ∫∑?≤ ∈
=
|)())((1)())((1|sup )()(
1
xdxWfcxdxWfcpc i
i
i
i
Lf
m
i
ii nm
∑
=
≤≤
m
i
ii cpc
1
),(),( nmrnmr,
这里我们用了,对于 Lf ∈,有 LWfc i
i
∈)(1 o,而 c 为 ic 们中的最大者,一旦证明了压缩性,那么 1?W 作为完备距离空间 ( M,r ) 的迭代映射的不动点就存在唯一,而且可以用迭代逼近,这个不动点就是我们所要的不变分布 函数 ),
遍历定理中关于压缩映射的条件可以减弱为 如下的 平均压缩性,这就是 Barnsley 与
Elton 在 1988 年得出 的 平均压缩条件下的遍历定理
定理９．３５ 若随机映射
)()1(
1
,,
,,
m
m
WW
pp
L
L 满足平均压缩条件,存在 1<c,使对于任
259
意 Xyx ∈,有
∏
=
≤
m
i
p
ii yxcdyWxWd
i
1
),())](),(([,(9,41)
那么,对于由独立同分布的如上的随机映射生成的随机的迭代函数系统，遍历定理的所有结论仍然成立,
(证 明 略 ),
(以上的不动点 (不变分布函数 ),作为图 像有某些自相似性,Barnsley 等人进一步希望用某种不动点来描述自相似的图象,为此提出了递归迭代函数系统 (Recursive Iterative Function System),他用一个 Markov
转移矩阵 P 来代替随机迭代函数系统中的随机向量 ),,( 1 mpp L,即记以 P 为转移矩阵的时齐的 Markov
链为 }{ nh,再 定义
)()(1 nn nW xx h=+,
此时 }{ nx 不再是 Markov 链,但是 },{ nn hx 却是 Markov 链,同样可以得到它的不变分布 函数,
另一方面,我们还可以在窗口上的全体概率测度的集合上,定义一个压缩映射得到不动点,而且它与
上面定义的是一样的,进一步还可以考虑有多个递归迭代函数的系统,得到这个系统的不动点 ),
4 统计中的 Bayes 方法与 图象的处理,分割 与重建
4．１ Bayes 统计要义
Bayes 方法的基本想法是,把概率函数 (包括分布密度和离散的分布概率 ) 中的 未知参数 q 当作随机变量 (或随机向量 ),
在 Bayes 方法中不再区分参数与随机变量．所以未知参数的分布的确认是最重要的．在抽取样本以前，就只能根据先验知识设置未知参数的分布，称为先验分布，或验前分布．在抽取样本以后，根据对于抽取到的样本的概率规律的了解，就可用 Beyes 公式把参数先验分布改进为后验分布，或验后分布，也称 Bayes 分布,
１．先验分布与后验分布 (Bayes 分布 )
设随机变量 x （或随机向量）的分布（或概率函数）为 ),( qxp （ q 可以是向量）．用
Bayes 的观点，把 q 看成随机变量后,),( qxp 就是：在 q 已知的条件下,x 的条件分布，
即
),()|(| qqqx xpxp =,
假定 q 的分布密度（或概率函数）为 )(Jf,称为 q 的 先验分布 (a priori distribution),,
于是 ),( qx 的联合分布为 ),()( qqf xp,在 J 取固定值的条件下,假若 nxx,,1 L 为采自
),( qxp 的独立随机样本,在 nxx,,1 L 已知的条件下,q 的条件分布记为
260
),,|( 1 nxxJj L,由 Bayes 公式可知它是 q 与 ( nxx,,1 L )的联合密度 ),()(
1
qxqf i
n
i
p∏
=
，除以 ( nxx,,1 L )的边缘密度，即
),,|( 1 nxxJj L
∫ ∏
∏
=
=
n
i
i
n
i
i
dp
p
1
1
),()(
),()(
JJxJf
JxJf
,(9,42)
称为 q 的 后验分布 ( a posteriori distribution),或 Bayes 分布,
2,已知损失函数是的的 Bayes 估计
假定对参数 J 的一个估计
^J
带来的损失可由一个 损失函数
^),( JJL
来度量．我们用
BayesE 表示在后验分布下取的数学期望,对于 J 的任意一个估计量
^J
,它在后验分布下的平均损失 ),(
^
JJLEBayes 称为 后验风险,使后验风险达到最小的估计,称为 Bayes 估计,记为
Bayes
^J
。即它满足
)],([min)],([ ^^ ^ JJJJ
J
LELE BayesBayesBayes =,( 9,43 )
3,Bayes 方法评述
Bayes 方法的长处是设置了一个先验分布，它可以发挥对未知参数的已有统计知识的认知．然而其难点也正在于如何选取合适的先验分布,这可以说是这种由从另一个角度认识 统计带来的一对矛盾,在实际应用中正是要巧妙地处理好这一对矛盾,
再则，由于在 Bayes 方法中，参数与随机变量是平等的，所以对于随机变量的估计也常用 Bayes 估计．特别 是 在数据不全的情形（缺失数据）的估计中，需要用测量数据来估计缺失数据，或者是在滤波问题中，或隐 Markov 模型中，需要用量测数据估计状态 （参见第
10 章与第 15 章）,
4．先验分布的参考取法
第一种取法出现在无任何先验信息的情形．具体地说
（ 1）在 J 取值于 ),( ∞?∞ 时，认为它是均匀的（这称为 Bayes 假设）,即先验密度为常数 (如果 q 可能取值区域面积为无限,则它不是真正的分布,而是所谓的,广义,分布，但是此时 (9,
43)仍然是一个分布，仍称为后验分布或 Bayes 分布 ),那么,Bayes 估计就变成为普通的最大似然估计,
（ 2）在 J 取值于 ),0[ ∞ 时，常用广义分布 J1 。
（ 3）在 J 取值于 ]1,0[ 时，常用广义分布 )1( 1 JJ? 。
第二种取法是用 最大 熵原则：取所谓的最大熵分布（参见第 10 章），它是在已知的限制条件下，信息掌握得最少的分布，说一句白话,就是“最吃亏的”分布,在对于参数完
261
全没有限制的情形下,广义的 最大熵分布恰恰就是均匀分布,所以第二种取法也是第一种取法的延伸。又若参数是一个概率向量（或者要估计的是一个取值于,概率向量组成的概率空间,的随机向量,这里
,n 维概率向量组成的概率空间” }1,1,,0:),,{( 111 =++≤≤=? nnn pppppp LLL ），
那么这时的最大熵分布是 Dirichlet 分布。这就是在生物信息论中常用 Dirichlet 分布作为先验分布的原因,
经验知识对于粗估先验分布常常很有用,例如,若未知的是英文字母或汉字,人们可以根据字母或汉字出现的经验频率得到先验分布的统计估计,
第三种取法是采用 Jeffreys 原则：取 先验密度为 )(det)( qqf IC=,其中 C 是常数，
)(qI 是信息矩阵，它的分量 )),,((),),(ln()( 1
2
k
ji
ij
xpEI qqq
qq L=
=,
第四种取法是统计中的所谓共轭分布方法,如果参数的先验分布与它关于分布 ),( qxp
的后验分布属于同一种类型，那么此种类型的分布就称为是 ),( qxp 的 共轭分布,人们认为用共轭分布作为先验分布有其合理性，因为这时先验分布与后验分布属于同一种类型,
常见分布的共轭分布有,Beta 分布是二项分布 的共轭分布,也是几何分布 的 共轭分布 ;
Gamma 分布是 Poisson 分布的共轭分布，也是指数分布的共轭分布；正态分布是已知方差的正态分布的共轭分布； 双侧截尾正态分布是已知方差的双侧截尾正态分布的共轭分布 ; 逆
Gauss 分布 是已知 均值 的正态分布 ( 包括 双侧截尾 情形 )的共轭分布 ; Perato 分布
（ )0,(),()( ),[)1( >?= ∞+? axIxaxp a gg ga ）是均匀分布的共轭分布 等,
在对参数没有任何先验知识时，在上面四种取法中，到底用那个更为合适？一般地说，
并没有确切的答案．这是 Bayes 方法的不足 之 处．实际上这时只能随 问题的不同作具体的考虑,
5,Bayes 估计再访
通过 Bayes 分布来估计参数或状态，有 3 种常见的方法,
（１）对给定的损失函数，求使 后验风险最小的估计,
（２）用后验分布的数学期望,即 后验期望估计,
（３）用后验分布达到最大值的点来估计，称为 最大后验估计,
与一般数理统计类似地还可以用后验分布作区间估计,
6,Bayes 策略
定义９,３６ 在统计抉择中，通常会涉及样本，参数与行动．参数 q 是被考虑成随机变量的,将 可能采取的行动全体组成的集合 记为 A．假定在对象的状态分布的参数为 q 时，
如果采取行动 )( Aa ∈,就会导致 损失 ),( aL q,而采取的动作 a 是依赖于被抽取到的样本
nxx,,1 L 的,我们将 它记为 ),,((),( 1 naa xxxx L==,这个“函数”（从样本到动作的对应） )(?a 就称为一个 策略 。 在样本为 x 时，策略 )(xa 带来的损失为 ))(,( xq aL,它 是样本的函数,
262
定义９,３７ 在 q 固定的时的平均损失为
xdxpxpxaLaELaR n ),(),())(,())(,())(,( 1 qqqxqq L∫==,(9,44)
称为 策略 )(?a 关于 q 的风险函数 。它对于先验分布 )(qf 的平均
=?))((ar ∫? qqfq daR )())(,(
(9,45)
称为 策略 )(?a 对于先验分布 f 的 Bayes 风险 。使先验分布 f 的 Bayes 风险达到最小的策略，
称为 对于先验分布 f 的 Bayes 策略 。
定义９,３８ 动作 )(xa 带来的损失 ))(,( xq aL 关于 q 的后验分布的期望，则称为 策略
)(xa 的后验风险 （就是 动作 )(xa 的后验风险），其表达式为
qxqxqxJ dpaLaLEBayes )|())(,())(,( ∫=,(9,46)
其中 )|( xqp 为参数 q 的 后验密度,即后验分布的密度,
这里需要分清楚,策略的风险函数是对状态变量作的平均,即对于一切可能的样本作平均，所以它定义在整个策略类上，而策略对先验分布的 Bayes 风险，是策略的风险函数关于先验分布的平均．而后验风险是，在给定样本下，策略对参数的后验分布作的平均．从数学处理上,在求后者时样本是固定的，所以要简单得多．这就是使用后验分布的好处，当然，后验分布的计算也常常并不简单,
需要注意，,Bayes”的冠名，通常用来修饰后验性质，例如,Bayes 分布就是后验分布．但是，在 Bayes 风险的定义 中，却被用来冠名先验性质,这似乎在冠名上有些混淆，
但是因为有下面的定理，所以并无实质上的不妥,
定理９．４０ (Bayes 策略基本定理 )
一个策略 )(*?a 是 Bayes 策略的充要条件是,对于任意样本 x,)(* xa 都使后验风险达到最小,
（此定理的证明需要较深的知识,从略）,
这个定理说明了,Bayes 策略在任意样本 x 上的取值，就是 x 的最小后验风险动作．于是最小 Bayes 风险在本质上仍反映了后验性质,
由此定理可见，只要算出了后验风险,再求得 Bayes 策略就并不困难。
Bayes 估计可以纳入 Bayes 策略的框架．再者，假设检验也可以纳入 Bayes 策略的框架，
因而 Bayes 统计方法的重点就归结为 Bayes 分布的计算。
在实际问题中,Bayes 分布的计算，常常要通过随机模拟得到,特别是 Bayes 分布中的分母的计算，常常涉及巨大的工作量，较为有效的数值计算方法，就是 Markov 链 Monte Carlo
方法（ MCMC）,
4,２ Bayes 方法在图像中的应用与 观测量不是状态变量时的参数估计
图象处理包括,分割 ( Segmentation),图象压缩,图象重建,图象分类 (如利用相似性对纹理图分类 ),图像识别,聚类等等,
263
图象在传输过程 先 需要进行压缩,然后由接收方 再将 压缩图象信息恢复,所以,图象处理主要为,图象压缩与图象恢复,而图象分割常用在,例如,自然场景的分割,在图象恢复,
图像分割的各种方法中,Bayes 方法是很重要的一个,它归结为在把接收到的信息作为已知的 条件下,对描述图象的一些参数的重估,
Bayes 聚类方法,就是根据观测 (从传输渠道得到 ),在轮廓图 (或灰度图 )全体组成的集合中,选取一个在此观测的 条件 下,Bayes 风险最小的轮廓图 (或灰度图 ),作为恢复的图象,
这里损失函数
^),( JJL
的取法是十分重要的,它必须取得使估计的结果较为稳健 (robust),
即更多地与分布的定性行为相关,而对具体分布的表示式的敏感性尽量地小,这里不能一味简单地搬用模式识别中的现成结果,而必须考虑到,标的图象,的来 源,结合它所在学科的特点,具体情形具体分析,
图像的先验分布 )(ap (Priori distribution),其中 a 为组态,一般可以写成 Gibbs 分布的形式 (只要它处处为正 ),
)()( aap He?=,
其中
)(ln)( apa?=H,
称为 能量函数 或 验前质量指标,
图像的信息特征中,一般还有代表边缘信息的 ‘边能量,( 对应于,边过程,),各部分的能量强度 的不同常表现为某些加权参数 jq 的大小,例如,图象恢复的模型可以取
∑
=
=
3
1
}|()|(
j
jjHH baqba,
这是一个参数建模问题,就是要确定参数,即要确定 H 为能量函数的 Gibbs 分布的参数 jq
们,在观测 )()1(,,mbb L （接收到的信息）条件下,作最大似然估计或 Bayes 估计,此时的联合分布 ( 即似然函数 ) 为
∏
∑=?
∑
∑
=
=
=m
k H
H
j
kj
j
k
j
e
eL
1 )|(
)|(
321 3
1
)(
3
1
)(
),,(
g
bgq
baq
qqq,
这里既有待估参数 q,又有未知组态 a,要同时估计它们,一种有效的方法是将在第 15章中介绍的 EM 算法,即用 EM 方法作迭代近似,对于格点数太多的情形,Geman 兄弟提出了一种简化的修正算法,可以大大地减少计算量,其基本思想为,一旦只要得到目标指向局部极大的
)(^ n
q 后,就按 ),|(
)(^
bq
n
P? 取样
)(^ n
a,然后按下面的方法更新
)1(^ +n
q,
任取
)1(^ +n
q 满足 (P
)(^ n
a ),|
)1(^
bq
+n
),|(
)(^)(^
bqa
nn
P>,
264
在实际估计中,对于参数的先验知识的利用是至关重要的,例如,参数所在的位置或区域,
Gibbs 场的 Bayes 分布纹理分析时,主要是根据纹理的特征,选取合 理的 Gibbs 分布,即选取能量函数,常用的能量函数有 Graffigne 的?Φ 模型,
∑?
+
=
相邻),(
2
2||
1
1)(
yx yx
xyH
d
aaqa,
在模型选取好以后,纹理图的聚类就化为：通过观测到的一系列标准样品,进行学习训练以确定 (估计 ) 参数 q,再用 Bayes 聚类或其它聚类法来决定此纹理图形属于那一类,
图像识别的融合原则
在图象识别时,有时可以利用不同的识别器的融合 (fusion),这时可以采取的原则,
例如有,多数原则,权威原则,民主原则 (以相似度 作 加权 )等,
习题 9
1,设 nw 独立同分布,且具有正连续密度 )( xr,又 bxga ≤≤< )(0,
)1(mod)(1 nnn wf +=+ xx,nnn wf +=+ )(1 hh,那么连续状态 Markov 链 nn hx,都是指数 1L 遍历的,
2,设 nw 独立同分布,且具有正连续密度 )( xr,假定 0>yf,且存在 0)( ≥yg,使在 x 固定时
),( yxf 对 y 的反函数,记为 ),(1 yxf?,满足 )(),(),(( 1 ygy yxfyxf≥?r,而 Markov 链 nx 满
足 ),(1 nnn wf xx =+,那么 nx 是指数 1L 遍历的 Markov 链,
3,设 nw 为独立同分布列，具有分布密度 r,证明
110011,,,),,(+?+ ==+= kknknnn xxwf xxxxx LL
是 k 阶 Markov 链,意即
),,,|( 0011 xxxAP nnnkn ===∈+ xxxx L )|( xAP nkn =∈= + xx,
