配分常数Z时的求和项太多,从而常会出现大量被求和的项的值小于计算机误差的精度, 以致使计算在实际上失去可能.用 Markov链 Monte carlo方法(MCMC)可以解决这个困难 即先近似地得到 Gibbs分布,再进行参数估计.所以MMC是一个重要的计算手段 MCMC方法就是构造以 Gibbs分布为可逆分布的 Markov链的转移矩阵.因为这里的随机 场的 Gibbs分布远比 Ising模型的 Gibbs分布复杂,我们只构造离散时间的 Markov转移矩 阵,而且还希望如 Glauber动力学那样,每次只在一个格点上改变(称为按 Gibbs方式转移) 最后还要保证由这样得到的转移矩阵的各行都收敛于 Gibbs分布 为了使构造更为透明,我们先对任给的格点子集J定义如下的“部分转移” B)={Z en(v)(a,B在/外坐标相同) (其它) 其中β,as表示:把组态a中在J处的坐标换成组态B中的相应坐标,而Z(a)是一个归 化常数(“局部的配分常数) ∑e (9.5) 这种“部分转移”表示组态间只在位置J外的坐标相同时才可能转移 在此基础上我们定义转移矩阵如下:记G中的全部格点的一个给定的排序为 {x1…,xa,(G|=N1N2).任给组态a,B∈S,定义 P,=(P,(a, B))a,Bes,P= 则P是一个转移矩阵.注意由于P中的转移只在一个格点x,上进行,由P构造的 Markov 链的一次转移,可以视为G|次连续转移的结果,其中每次转移只在一个格点上进行,而且 要做遍所有的格点 下面的遍历收敛定理是 Markov链的 Monte carlo方法的理论基础.而 Markov链 Monte Carlo方法主要用于获得 Gibbs分布的样本,以便进一步得到Gibs场的各种统计平均量. 对于任意格点x∈G,任意组态a∈S,记 Hu(r,aGuri)=m,=min n, Hu(nraGyr) 那么 (Hu(,:)-m,I Pir(a, B -Hu(B! -m, e max( Hu(B)Hu(y)BGx=yGx)A G 其中δx为能量函数HU在格点x上的振幅,Δ为Hu的振幅,其确切定义为
234 配分常数 Z 时的求和项太多，从而常会出现大量被求和的项的值小于计算机误差的精度， 以致使计算在实际上失去可能．用 Markov 链 Monte Carlo 方法 (MCMC)可以解决这个困难， 即先近似地得到 Gibbs 分布, 再进行参数估计. 所以 MCMC 是一个重要的计算手段. MCMC 方法就是构造以 Gibbs 分布为可逆分布的 Markov 链的转移矩阵. 因为这里的随机 场的 Gibbs 分布远比 Ising 模型的 Gibbs 分布复杂, 我们只构造离散时间的 Markov 转移矩 阵, 而且还希望如Glauber动力学那样, 每次只在一个格点上改变（称为按 Gibbs方式转移）, 最后还要保证由这样得到的转移矩阵的各行都收敛于 Gibbs 分布. 为了使构造更为透明, 我们先对任给的格点子集 J 定义 如下的“部分转移” ï î ï í ì = - 0 ( ) ( , ) ( ) 1 ( , ) ( ) \ 其它 e 在J外坐标相同 p Z HU J S J J J a b a b a b a , (9. 4) 其中 b JaS \ J 表示: 把组态a 中在 J 处的坐标换成组态 b 中的相应坐标, 而 (a) ZJ 是一个归 一化常数(“局部的”配分常数): å - = J HU J S J J Z e b b a a ( ) \ ( ) . (9. 5) 这种 “部分转移”表示组态间只在位置 J 外的坐标相同时才可能转移. 在此基础上我们定义转移矩阵如下：记 G 中的全部格点的一个给定的排序为 { , , ),(| | ) 1 | | G N1N2 x x L G = . 任给组态a,b Î S , 定义 P J J S p = a b a ,bÎ ( ( , )) , P = Õ= | | 1 G k P ( xk } , (9. 6) 则 P 是一个转移矩阵. 注意由于 P ( xk } 中的转移只在一个格点 i x 上进行, 由 P 构造的Markov 链的一次转移,可以视为| G |次连续转移的结果, 其中每次转移只在一个格点上进行, 而且 要做遍所有的格点. 下面的遍历收敛定理是 Markov 链的 Monte Carlo 方法的理论基础. 而 Markov 链 Monte Carlo 方法主要用于获得 Gibbs 分布的样本, 以便进一步得到 Gibbs 场的各种统计平均量. 对于任意格点 x ÎG ,任意组态a Î S , 记 ( ) min ( } HU xa G\{ x} mx x HU hxaG\{ x} g h D D = = . (9. 7) 那么 ( , ) p{x} a b åÎ - - - - = y G H m H m U y G y x U x G x x e e [ ( ) ] [ ( ) ] \{ } \{ } b a b a - -D D - - = ³ = ³ e G e G G e x HU HU G x G x | | 1 | | 1 | | max( ( ) ( ): ) \{ } \{ } d b g b g , 其中 x d 为能量函数 HU 在格点 x 上的振幅, D 为 HU 的振幅, 其确切定义为
δ,= max a(x2=a I Hu(a)-Hu(B) 由此推出 C(Px)≤l-G| min aB Pirl(a,B)≤1 从而得到 (P)=C(Px1)…C(Px)≤(1 (9.10) 定理9.12(遍历收敛定理)以上定义的转移矩阵P的 Dobrushin收缩系数满 C(P)≤1-e(△为H的振幅) 记以P为转移矩阵的 Markov链为n,则它以Gibs分布兀为可逆分布,而且有 证明显见r(a)p(a,B)=r(B)P(B,a)·由P=(p(,B)a,Bes的定义立得 丌(a)p(a,B)=(B)p(B,a).因此兀为可逆分布.又由C(P)<1,用 Dobrushin定理可 知极限成立 [注]G-格点的排序{x,…,xa)并不是实质的我们可以用一个G上取值全正的预选概率分 布g(x),(即g(x)>0∑g(x)=1例如均匀分布)的随机数作随机选取来代替.这时转移矩阵P应 作相应的改动,代替P1x(a,B)的是,以g(x)的概率对只在x坐标不同的组态间进行转移,即令 =(P),P=(p(ax,B) 其中 ,B)=8(x)(x,B)(若o=Bm:) (其它情形) 这时我们仍有 1 1.5 Gibbs分布的模拟退火 本段讨论求得达到组态空间上能量函数最小值的组态的概率方法,就是用能量函数的 Gibs分布作模拟退火的方法.考虑依赖于温度Tn的Gibs分布,以它们作为不变分布构造 各步为非时齐Gibs转移矩阵,由此确定的一个非时齐 Markov链{n}.我们可以用它作模 拟退火,使得当n充分大以后,此非时齐的 Markov链n以大概率落在能量函数最小的组态 5
235 max | ( ) ( )| \{ } \{ } d x a b HU a HU b G x G x = = - , x G x d Î D D =max . (9. 8) 由此推出 -D C P £ - G p £ - e ( {x} ) 1 | | min a ,b {x} (a,b) 1 . (9. 9) 从而得到 | | { } { } ( ) ( ) ( ) (1 ) 1 G x x C P C P C P e M -D = L £ - | | 1 G e -D £ - . (9. 10) 定理 9. １２ （遍历收敛定理） 以上定义的转移矩阵P 的 Dobrushin 收缩系数满 足 | | ( ) 1 G C P e -D £ - ( D 为 HU 的振幅). 记以 P 为转移矩阵的 Markov 链为 n x , 则它以 Gibbs 分布p 为可逆分布, 而且有 P p n n T ¾ ¾®1 ®¥ . 证明 显见 p (a) (a,b ) p (b ) (b,a) pJ = pJ . 由 P S p = a b a,bÎ ( ( , )) 的定义立得 p (a) p(a, b ) = p (b ) p(b ,a) . 因此p 为可逆分布. 又由 C(P) < 1, 用 Dobrushin 定理可 知极限成立. [注] G -格点的排序 { , , ) 1 |G| x L x 并不是实质的. 我们可以用一个 G 上取值全正的预选概率分 布 g (x) , (即 ( ) > 0,å ( ) = 1 xÎG g x g x 例如均匀分布) 的随机数作随机选取来代替. 这时转移矩阵 P 应 作相应的改动，代替 ( , ) { } a b xk p 的是, 以 g (x) 的概率对只在 x 坐标不同的组态间进行转移, 即令 P ( D = ~ P | | ) G , ~ P = p a b a,bÎS ~ ( ( , )) , 其中 î í ì = = 0 ( ) ( ) ( , ) ( ) ( , ) { } \{ } \{ } ~ 其它情形 p x 若 G x G x g x p a b a b a b . 这时我们仍有 P p n n T ¾ ¾®1 ®¥ . １. ５ Gibbs 分布的模拟退火 本段讨论求得达到组态空间上能量函数最小值的组态的概率方法, 就是用能量函数的 Gibbs 分布作模拟退火的方法. 考虑依赖于温度Tn的 Gibbs 分布, 以它们作为不变分布构造 各步为非时齐 Gibbs 转移矩阵，由此确定的一个非时齐 Markov 链{ }n x . 我们可以用它作模 拟退火, 使得当n 充分大以后，此非时齐的 Markov 链 n x 以大概率落在能量函数最小的组态
所组成的集合中 更具体地说,对于如下的依赖于温度Tn的O一邻位势Gibs分布 -Hula) -Hu(a) 相应地定义 pj(a, B)=z.(a) (a,B在/外坐标相同) (9.11) (其它) -Hu(Busy) (a)= ∑ 再令 1(=(P1(a,B)s,P=∏P,(9.12) TEG 则P(是一个转移矩阵.我们将以它为非时齐转移矩阵列所确定的非时齐 Markov链记为 {n} 定理9.13( Gibbs分布的模拟退火的收敛性定理) 如果温度T满足如下的 Geman- Geman条件:对0<T√0,存在n0使n≥n0后有 Tn≥ (9.13) In 其中Δ是能量函数Hu(a)的振幅 A=maxes Hu(a)-m H() 那么在任意初始分布下,以P为第n步转移矩阵列的非时齐的 Markov链{n},在时 刻n的绝对分布n有极限 n→ 其中丌)为集合{a∈S:H(x)=min}上的离散均匀分布.从而还有 lm n P(Hy(sn)=minges hy(a))=1 证明我们验证 Dobrushin- Isaacson- Madsen定理的条件成立.条件(A.1)的验证与第 8章3.2段的模拟退火定理中的证明完全一样.下面验证 Dobrushin条件(A.2)’.这里 定理9.12中的不等式变为 36
236 所组成的集合中. 更具体地说, 对于如下的依赖于温度Tn的¶ -邻位势 Gibbs 分布 ( ) 1 ( ) 1 ( ) a p a U n H T n n e Z - = , åÎ - = S H T n U n Z e a (a ) 1 , 相应地定义 ï î ï í ì = - 0 ( ) ( , ) ( ) 1 ( , ) ( ) 1 , ( ) \ 其它 e 在J外坐标相同 p Z U J S J n H T n J n J a b a b a b a , (9. 11) å - = J U J S J n H T n J Z e b b a a ( ) 1 , \ ( ) . 再令 P ( ( ) { } D = n x S n = p x a b a,bÎ ( ) { } ( ( , )) , P = (n) ÕxÎG P ( ) { } n x , (9. 12) 则 P (n) 是一个转移矩阵. 我们将以它为非时齐转移矩阵列所确定的非时齐 Markov 链记为 { }n x . 定理９.１３（Gibbs 分布的模拟退火的收敛性定理） 如果温度Tn满足如下的 Geman-Geman 条件: 对0 < ¯ 0 Tn , 存在n0 使n ³ n0后有 n G Tn ln | | D ³ , (9. 13) 其中D 是能量函数 (a) HU 的振幅: D max (a) min (a) a a D = ÎS HU - ÎS H , 那么在任意初始分布 m0下, 以P (n) 为第n 步转移矩阵列的非时齐的 Markov 链{ }n x , 在时 刻n 的绝对分布 mn 有极限 ® (¥) mn p , 其中 (¥ ) p 为集合{a Î : (a) = min} S HU 上的离散均匀分布. 从而还有 lim n®¥ P(HU (xn ) = min aÎS HU (a)) = 1. 证明 我们验证 Dobrushin-Isaacson-Madsen 定理的条件成立. 条件(A.1)的验证与第 8 章 3. 2 段的模拟退火定理中的证明完全一样. 下面验证 Dobrushin 条件(A. 2)’. 这里 定理９.１２中的不等式变为
其中△为能量函数HU的振幅.于是对于PP+1)…P,我们有 C(PP)…P)≤C(PC(P+)…C(P)≤∏(1-e) 由无穷乘积的理论,当n→∞时,上面右方趋于0的充要条件为级数∑e是发散的 但是由 Geam- Geam条件,我们有 所以 Dobrushin条件(A.2)成立 [注1] Gibbs场的模拟退火在图象的清污(图象的滤波)中的应用 设图象为组态a,而在传输过程这受到了污染,使得接受方观测到的图象为组态B.需要用数学方法 进行滤波,以便恢复图象的清晰度,我们把在观测为组态B的条件下,状态a的分布,记为P(a|B) (就是后验分布 Posteriori distribution,即 Baves分布).一般它也可用写成Gbbs分布的形式 P(aB) Z(B) 其中H(a|B)称为条件能量函数.典型的例子如 H(aB)=0∑(a(x)-a(y)2+∑a(x)-B(x)2 其中第一项强调了图形的光滑性,第二项为保真度.滤波的目的就是要找组态α,使条件能量函数达到最 H(a|β)= min h(a|B 这个最佳组态O可以通过模拟退火方法求得 又例如,在边缘检测中,β就是离散化点上的象素,而α就是“边随机链”的一个样值,与滤波问 的不同的是,此时的后验图形只是先验图形的一部分,即边的位置.在复杂的图形处理技术中,例如,X断 层摄影术的处理,情况就远为复杂.适当的利用分割技术可以得到较好的恢复效果 [注2]图象的处理有许多数学方法, Gibbs分布方法只是其中的一种.另一种更为常用的方法是,把 图象的一些特征参数作为观察量(这同时也达到了图象信息压缩的目的.用它们来估计 Markov链的状态 就是重建图象.这就是隐 Markov方法(M),一般地说,如果特征参数取得合适,HM方法是十分有效的. 例如,纹理图象分析方法( Texture Analysis)提取特征参数.其它的方法还有:下节将介绍的随机函数迭
237 - D £ - Tn G n C P e | | ( ) ( ) 1 , 其中 D 为能量函数 HU 的振幅. 于是对于 (j) (j 1) (n) P P LP + , 我们有 £ + ( ) ( j) ( j 1) (n ) C P P LP £ + ( ) ( ) ( ) ( j) ( j 1) (n) C P C P LC P (1 ) | | | | n T G G k j e D - = Õ - . 由无穷乘积的理论, 当n ® ¥时, 上面右方趋于 0 的充要条件为级数å D - n T G n e | | 是发散的. 但是由 Gemam-Gemam 条件, 我们有 å D - n T G n e | | å - ³ n n e ln = å = ¥ n n 1 . 所以 Dobrushin 条件(A. 2)’成立. ［注 1］ Gibbs 场的模拟退火在图象的清污(图象的滤波)中的应用 设图象为组态a , 而在传输过程这受到了污染, 使得接受方观测到的图象为组态 b . 需要用数学方法 进行滤波, 以便恢复图象的清晰度. 我们把在观测为组态 b 的条件下, 状态a 的分布, 记为P(a | b ) (就是后验分布 Posteriori distribution, 即 Bayes 分布). 一般它也可用写成 Gibbs 分布的形式: ( | ) ( ) 1 ( | ) a b b a b H e Z P - = , (9. 14) 其中 H (a | b )称为条件能量函数. 典型的例子如 H (a | b ) å å ζ = - + - ( ) 2 2 2 ( ( ) ( )) 2 1 ( ( ) ( )) x y x x y a x b x s q a a , 其中第一项强调了图形的光滑性, 第二项为保真度. 滤波的目的就是要找组态 ^ a , 使条件能量函数达到最 小: ( | ) min ( | ) ^ H a b = a H a b . 这个最佳组态 ^ a 可以通过模拟退火方法求得. 又例如, 在边缘检测中, b 就是离散化点上的象素, 而a 就是 “边随机链”的一个样值, 与滤波问题 的不同的是, 此时的后验图形只是先验图形的一部分, 即边的位置. 在复杂的图形处理技术中, 例如, X 断 层摄影术的处理, 情况就远为复杂. 适当的利用分割技术可以得到较好的恢复效果． [注 2] 图象的处理有许多数学方法, Gibbs 分布方法只是其中的一种. 另一种更为常用的方法是, 把 图象的一些特征参数作为观察量(这同时也达到了图象信息压缩的目的. 用它们来估计 Markov 链的状态 , 就是重建图象. 这就是隐 Markov 方法(HMM), 一般地说, 如果特征参数取得合适, HMM 方法是十分有效的. 例如, 纹理图象分析方法(Texture Analysis)提取特征参数. 其它的方法还有: 下节将介绍的随机函数迭
代系统方法.除此以外,还有小波分析等非概率方法, [注3]当格点集G为无限集时,组态的数目就不再是可数的了.此时就应归入连续状态的 Markov 场的范畴.在统计物理中应用的模型,就是这种模型.有时自旋空间S也可以是连续的,例如,圆周或其 它集合.这时就可能存在多个不变分布,称为有相变,因为每个分布代表了统计物理中质点的一个宏观 相”.统计物理中也正是对这种相变的存在性更感兴趣. [注4]若每次发展的修改的格点集合J不是单点集,且各次的修改集可以不同,那么,此时的 Markov 链未必有遍历定理,这时并不能保证此 Markov链5n的时间平均收敛到不变分布,因此,即使n很大,也 不能认为5n的样本近似为不变分布的样本 2时间离散状态连续的 Markov链 概率空间再访 设g2={0:O-基本事件},7是样本空间的某些子集组成的集合类(事件体),满足 对于任意AAn∈?(n21),恒有UAn∩A1,A=92-A∈只 (这样的事件体?代表了可以通过理论推算,或实际测量能够得到概率的全体随机事件).在 概率理论中?称为Ω上的一个σ-代数.(Ω,7)称为概率空间.于是概率P就是定义在 σ-代数?上的一个非负函数 定义9.14( Borel集与 Borel函数) 设是R的某些子集组成的一个-代数,且满足以下条件 (1)R4中的任意开集G∈ (2)任意一个有R的某些子集组成的σ-代数另,只要它含有一切开集,则就有 那么因称为R“上 Borel集合类,其中的集合称为Bore|集.即 Borel集合类是以开集为元 素的最小σ-代数.(以前在本书中的所谓“常见的”集合,均指 Borel集).Rl上的 Borel 集称为一维 Borel集 再设w为R上的一个满足下述条件的函数类 1)对极限封闭:若∫n∈w且fn(x)→f(x),则有∫∈ (2)任意连续函数∫∈w, (3)若R上一个函数类∠只要对极限封闭,而且包含全体连续函数,则有wC∠
238 代系统方法. 除此以外, 还有小波分析等非概率方法. [注 3] 当格点集 G 为无限集时, 组态的数目就不再是可数的了. 此时就应归入连续状态的 Markov 场的范畴. 在统计物理中应用的模型, 就是这种模型. 有时自旋空间 S 也可以是连续的, 例如, 圆周或其 它集合. 这时就可能存在多个不变分布, 称为有相变, 因为每个分布代表了统计物理中质点的一个宏观 “相”. 统计物理中也正是对这种相变的存在性更感兴趣. [注 4] 若每次发展的修改的格点集合 J 不是单点集, 且各次的修改集可以不同, 那么, 此时的 Markov 链未必有遍历定理, 这时并不能保证此 Markov 链 n x 的时间平均收敛到不变分布, 因此, 即使n 很大, 也 不能认为 n x 的样本近似为不变分布的样本. 2 时间离散状态连续的 Markov 链 2. 1 概率空间再访 设W ={w :w - 基本事件}，F 是样本空间W 的某些子集组成的集合类(事件体), 满足 对于任意 A, An ÎF (n ³ 1)，恒有 U I ¥ = ¥ = = - Î 1 1 , , n n n An A A W A D F. (这样的事件体F 代表了可以通过理论推算, 或实际测量能够得到概率的全体随机事件). 在 概率理论中 F 称为 W 上的一个s - 代数. (W , F ) 称为概率空间. 于是概率P 就是定义在 s - 代数 F 上的一个非负函数. 定义９.１４（Borel 集与 Borel 函数） 设 B 是 d R 的某些子集组成的一个s - 代数，且满足以下条件: (1) d R 中的任意开集GÎB , (2) 任意一个有 d R 的某些子集组成的s - 代数 G, 只要它含有一切开集, 则就有 B Ì G． 那么 B 称为 d R 上 Borel 集合类, 其中的集合称为 Borel 集. 即 Borel 集合类是以开集为元 素的最小s - 代数. (以前在本书中的所谓 “常见的”集合, 均指 Borel 集). 1 R 上的 Borel 集称为一维 Borel 集． 再设 M 为 d R 上的一个满足下述条件的函数类: (1) 对极限封闭: 若 f n ÎM 且 f (x) f (x) n ® , 则有 f ÎM, (2) 任意连续函数 f ÎM, (3) 若 d R 上一个函数类 L, 只要对极限封闭, 而且包含全体连续函数, 则有 MÌ L．