龚光鲁,钱敏平著应用随机过程教程一与在算法和智能计算中的应用 清华大学出版社,2003 第15章与数据建模有关的几个算法 1EM算法一隐状态变量分布中参数的最大似然估计 EM算法的基本想法 在数据资料不全时,由已有的资料Y估计缺失变量X,或在观测到的资料Y并不是状 态变量时,估计状态变量X时(或者估计其分布密度f(9,x)中的未知参数9),就与古典统 计很不相同,此时需要用观测到的资料Y,同时估计X与未知参数 这样的估计将面临如下困难:如果把在参数9下的期望记为E,那么,在估计状态 变量X时,估值当然应该用条件期望X=E(X|Y)(如果在Y=y及9的条件下,X的 Bas分布密度∫(x1,9)为已知,则也常用 Bayes估计x=x(x1,9)x).然而这 时就需要知道参数9的值;另一方面,为了知道9,又必须先知道X的估值X(作为状态 的已知样本值).这样,估计状态与估计未知参数之间是耦合的 在统计中通常对付这类困难的解耦方法是:假定一个已知,迭代地分别交替估计它们中 的另一个,直至稳定.此类算法通称为酬M算法,较为确切的表达是: (1)设置初值9 (2)(E-步骤)对n≥0,令X(m=E(X|Y)(或用 Bayes估计 f(x r, 9, dx ) (3)(M-步骤)(修正⑨的估计)取9n+1使之满足: log f(9m I, X )=Maxg log f(9, x) 其中E-步骤为取条件期望( Expectation),而M-步骤为取最大( Maximu 这种交替 迭代的方法,称为简单的EM方法 这个算法的构思很简单,但计算量过大,且一般很难看出是否稳定.为了克服这个缺 点, Dempster, Laird和 Rubin提出了直接递推估计9的想法(仍旧称为EM方法),这种 经过本质改进后的方法,至少在直观上看起来有稳定趋势 1.2 Rubin算法 假定(x,Y)具有联合分布密度f(9,x,y) Rubin算法的核心构思为:直接使用状态变量Y的分布密度g(9,y)代替(X,Y)的分 419
419 龚光鲁, 钱敏平著 应用随机过程教程 – 与在算法和智能计算中的应用 清华大学出版社, 2003 第 15 章 与数据建模有关的几个算法 1 EM 算法 – 隐状态变量分布中参数的最大似然估计 1. 1 EM 算法的基本想法 在数据资料不全时，由已有的资料Y 估计缺失变量 X ，或在观测到的资料Y 并不是状 态变量时,估计状态变量 X 时（或者估计其分布密度 f (J, x) 中的未知参数J ），就与古典统 计很不相同，此时需要用观测到的资料Y ，同时估计 X 与未知参数J . 这样的估计将面临如下困难: 如果把在参数J 下的期望记为 EJ , 那么, 在估计状态 变量 X 时，估值当然应该用条件期望 ( | ) ^ X = EJ X Y ( 如果在Y = y及J 的条件下, X 的 Bayes 分布密度 f (x | y,J) 为已知, 则也常用 Bayes 估计 ò X = x f (x | Y, )dx ^ J ). 然而这 时就需要知道参数J 的值; 另一方面，为了知道J ，又必须先知道 X 的估值 ^ X (作为状态 的已知样本值). 这样, 估计状态与估计未知参数之间是耦合的. 在统计中通常对付这类困难的解耦方法是: 假定一个已知，迭代地分别交替估计它们中 的另一个, 直至稳定．此类算法通称为 EM 算法, 较为确切的表达是: (1) 设置初值 J0； (2) (E-步骤) 对n ³ 0, 令 ( | ) ^ ( ) X E X Y n n = J (或用 Bayes 估计 ò X = x f x Y dx n n ( | , ) ^ ( ) J )； (3)（M -步骤）（修正J 的估计） 取Jn +1使之满足： log ( , ) log ( , ) ^ ( ) ^ ( ) 1 n n f Jn+ X = MaxJ f J X , 其中 E-步骤为取条件期望 ( Expectation ), 而 M-步骤为取最大（ Maximum ）. 这种交替 迭代的方法, 称为简单的 EM 方法. 这个算法的构思很简单, 但计算量过大, 且一般很难看出是否稳定. 为了克服这个缺 点, Dempster, Laird 和 Rubin 提出了直接递推估计J 的想法(仍旧称为 EM 方法), 这种 经过本质改进后的方法, 至少在直观上看起来有稳定趋势. 1. 2 Rubin 算法 假定(X,Y) 具有联合分布密度 f (J, x, y) Rubin 算法的核心构思为： 直接使用状态变量Y 的分布密度 g(J, y) 代替(X,Y) 的分
布密度∫(9,x,y),来求9关于g(9,y)的最大似然估计9.也就是求使logg(9,Y)达到 最大的9 由于g(9,y)在实际上并不好求, Dempster- Laird-Rbin利用了下述等式 f(9,x,y)=f(9,x|y)g(9,y) (15.1) 于是 log g(o, Y)=log f(, X, y)-logf(, XY) (15.1) 其右方是状态变量的对数似然密度与对数条件似然密度的差.对等式两边取关于Y的条件 期望得到 L(o)=log g(o, r)(=Eg(og g(o, r)r) g( Llogf(o, X, Y)IY-Es log [f(,X Q(q|9)-H(q|9). (15.2) 在观测资料是y的时候即Y=y的时候),Q(q|9)的表达式为 g(o|9)= log /(p, x,D)xr(0,xly)d 而H(|9)的表达式为 H(o19)= log/(o, x,yf(e, xly) =∫lgf1(X1)所(,x1yx 为了求L(9)的极大值点,我们将沿用第10章中对得到隐 Markov模型的模型参数估 计的 Baum-Welsh算法的思想.为此注意 L(q)-L(9)=[Q(q|9)-Q(9|9+[H(9|9)-H(q|9 (9,x|y) Q(o|9)-Q(9|9)+∫pog ≥[Q(q|9)-Q(9|9, (15.3) 其中用到第2项是相对熵的非负性.由(15.3)可知 只要Q(q|9)-Q(9|9≥0就有L(q)-L(9)≥0 420
420 布密度 f (J, x, y) , 来求J 关于 g(J, y) 的最大似然估计 ^ J ． 也就是求使log g(J,Y) 达到 最大的 ^ J . 由于 g(J, y) 在实际上并不好求, Dempster– Laird- Rubin 利用了下述等式 f (J, x, y) = f (J, x | y) g(J, y) ， (15. 1) 于是 log g(j,Y) = log f (j, X,Y) - log f (j, X |Y) ， (15. 1)’ 其右方是状态变量的对数似然密度与对数条件似然密度的差. 对等式两边取关于Y 的条件 期望得到 L(j) log g(j,Y) D = （ E (log g(j,Y) | Y) = J ） E log f (j, X,Y) | Y E log f (j, X | Y)] = J（［ ］ ）- J ［ =Q(j |J) - H (j |J) D . (15. 2) 在观测资料是 y 的时候(即Y = y 的时候), Q(j |J) 的表达式为 Q(j |J) f x y f x y d x X Y log ( , , ) ( , | ) | j q ò = , 而 H (j |J) 的表达式为 H (j |J) = f x y f x y d x X Y log ( , , ) ( , | ) | j q ò f x y f x y d x X Y X Y log ( , | ) ( , | ) | | j q ò = . 为了求 L(J) 的极大值点, 我们将沿用第 10 章中对得到隐 Markov 模型的模型参数估 计的 Baum-Welsh 算法的思想. 为此注意 L(j) - L(J) = [Q(j |J) - Q(J |J)] + [H(J |J) - H (j |J)] f x Y dx f x Y f x Y Q Q X Y X Y X Y ] ( , | ) ( , | ) ( , | ) [ ( | ) ( | )] [log | | | J j J j J J J ò = - + ³ [Q(j |J) - Q(J |J)] , (15. 3) 其中用到第２项是相对熵的非负性. 由(15. 3)可知: 只要 Q(j |J) - Q(J | J) ³ 0 就有 L(j) - L(J) ³ 0
所以,为了将9的较粗估计9,修改为较精的估计⑨n1,只需找9n1使 L(9n1)-L(9n)≥0.也就是只需要找n1使 Q(9n1|9n)-Q(9n|9n)≥0 (15.4) 于是我们就得到如下的 Dempster- Laird- Rubin的酬M算法,简称 Rub in算法或EM算法 (1)设置初值9 (2)(E步骤)对于n20,计算Q(|9,)=「g/(,xy)/(0,x1y)dx (3)(M-步骤)(修正9的估计)取⑨n+!使 Q(9m+1 9n)=Max @( 9n) (15.5) 由于将9n修正为9n时,Q(9n1|9)比Q(9n|9n)大,所以L(⑨n)是不减的,这就说明 9有收敛的趋势. Rubin证明了,在一定的条件下,9确实按概率收敛于某个9.于是我 们可以合理地把9作为分布的参数9的较佳的估计 注1L(9)是多峰函数的时候,9有可能是L(9)的局部最大值,甚至是鞍点,为了 纠正这种不足,一般可以从多个初始值出发,找到多个9值进行比较 注2所谓“缺失资料”可见之于两种情形:一种是观测资料的丢失,另一种是人为 地设置一些“后台操作的”辅助随机变量,称之为潜变量,将他们视为缺失的资料,即将他 们视为不能直接测量到的状态随机变量.这种潜变量的取法可以非常灵活,依赖于人们对于 问题的认识,经验的积累与对于技术掌握的成熟程度.从数学处理的角度考虑,最好能选取 潜变量X,使它与观测随机变量γ的联合分布是指数型分布(参见第1章),这时E-步骤就 很容易计算. 注3在离散随机变量的情形,E一步骤中的积分应该相应的改为求和。 以上算法也同样将隐马氏模型中的Baum- Welsh算法,纳入了EM算法的框架. 1.3EM算法的变通一广义EM算法 在实际计算中,又因为E-步骤的计算量很大,人们常常并不真去计算条件期望,而 是采用随机模拟.即:用在(9,})已知的条件下,用随机模拟得到X的 Bayes分布的若 干独立随机数,代入log∫(φ,Ⅺ)后作样本平均,用来代替Q(φ|9n)·这就大大地减少了 计算量.实践表明这样的简化,常可以得到相当满意的效果 同样,也因为M步骤的计算量也很大,人们也常常并不去算最大,而是任意找一个 φ,只要满足Q(φ|9)-Q(9|9n)≥0,就取φ为⑨n1·这样简化了的算法,称为广义 421
421 所 以 , 为 了 将 J 的较粗估计 Jn , 修改为较精的估计 Jn+1 , 只需找 Jn+1 使 L(Jn+1 ) - L(Jn ) ³ 0. 也就是只需要找Jn+1 使 Q(Jn+1 |Jn ) -Q(Jn | Jn ) ³ 0 . (15. 4) 于是我们就得到如下的 Dempster– Laird- Rubin 的 EM 算法, 简称 Rubin 算法或 EM 算法: (1) 设置初值 J0； （2）（E-步骤）对于n ³ 0, 计算 ( | ) Q j Jn f x y f x y d x X Y log ( , , ) ( , | ) | j q ò = ； （3）（M -步骤）（修正J 的估计） 取Jn +1 使 ( | ) ( | ) Q Jn+1 Jn = Maxj Q j Jn . (15. 5) 由于将Jn修正为Jn+1时, ( | ) Q Jn+1 Jn 比 ( | ) Q Jn Jn 大, 所以 ( ) L Jn 是不减的, 这就说明 Jn有收敛的趋势. Rubin 证明了，在一定的条件下，Jn确实按概率收敛于某个 ^ J ．于是我 们可以合理地把 ^ J 作为分布的参数J 的较佳的估计 注 1 L(J) 是多峰函数的时候， ^ J 有可能是 L(J) 的局部最大值，甚至是鞍点，为了 纠正这种不足，一般可以从多个初始值出发，找到多个 ^ J 值进行比较． 注 2 所谓“ 缺失资料”可见之于两种情形：一种是观测资料的丢失，另一种是人为 地设置一些“后台操作的”辅助随机变量，称之为潜变量，将他们视为缺失的资料，即将他 们视为不能直接测量到的状态随机变量. 这种潜变量的取法可以非常灵活, 依赖于人们对于 问题的认识, 经验的积累与对于技术掌握的成熟程度. 从数学处理的角度考虑, 最好能选取 潜变量 X ,使它与观测随机变量Y 的联合分布是指数型分布 (参见第 1 章), 这时 E-步骤就 很容易计算. 注 3 在离散随机变量的情形，E－步骤中的积分应该相应的改为求和。 以上算法也同样将隐马氏模型中的 Baum-Welsh 算法, 纳入了 EM 算法的框架． １．３ EM 算法的变通 – 广义 EM 算法 在实际计算中, 又因为 E-步骤的计算量很大, 人们常常并不真去计算条件期望, 而 是采用随机模拟. 即: 用在( ,Y) Jn 已知的条件下, 用随机模拟得到 X 的 Bayes 分布的若 干独立随机数, 代入 log f (j, X) 后作样本平均, 用来代替 ( | ) Q j Jn . 这就大大地减少了 计算量. 实践表明这样的简化, 常可以得到相当满意的效果. 同样, 也因为 M-步骤的计算量也很大, 人们也常常并不去算最大, 而是任意找一个 j , 只要满足Q(j |Jn ) - Q(Jn |Jn ) ³ 0, 就取j 为 Jn+1 . 这样简化了的算法, 称为广义
EM算法,常称为GEM算法 注]在应用EM方法时,为了避免M-步骤中求最大值的复杂计算,还可以采取其它的灵活的替代 方法.例如有 ECM算法〔 conditional max imum,cM),在多个参数的时候,例如9=(91,92)的情形,如下的交 替地求条件极值的最大化的方法,也常被用来代替M-步骤: (1)先取92满足 Q(9),92|9)=maxQ(91",92)9") (2)再取91满足 Q(9),92m)1(91",92")= maxi. o(9',92")1(9",92") 称之为ECM算法。一般地,ECM算法比M算法达到稳定的时间长 另有一种ECME算法(混合算法),它就是交替地使用FCM算法和M算法 2.在数据不完全时,用增补潜在数据,对参数的 Bayes分布作估计- Tanner-Wong的潜变量法 基本想法一估计后验分布 在数据资料不全(缺失数据),或观测到的资料Y并不是状态变量时,估计状态变量X的分布密度 f(9,x)中的未知参数9,与古典统计不同处是只能用观测到的资料Y,这时可以通过后验分布密度 p(9|Y)得到参数9的 Bayes估计.然而,后验分布p(9|Y)一般并不知道,就需要用观测到的资料y 对后验分布p(9|y)进行估计,这就是本段的目的 Tanner-Wong的想法是:在某些情形下,由观测数据Y可以通过条件分布取样的机制,来构造某 种”增补数据”(记为Z,称为潜变量, L atent Variable)的样本值.于是Z的这些样本值取自条件分布密 度P(=|y,9).P(=|y)称为预测分布.而在Y=y已知的条件下,潜变量的各个样本值彼此是条件 独立的 潜变量Z是观测不到的,如何选好它,最为重要.最简单的情形是:观测数据Y=Y1,2…,YN 为N个时刻的历史资料,其中Y1是一个m维向量,如果它的第1个分量是缺失的,而其它m-1维就 是状态变量.这时候潜变量就可取为那个缺失的一维变量 然而,一般情形远非如此简单.潜变量的选取是关系到能否有效地计算的关键, Tanner-wong选取 z的原则是,同时满足以下两个条件(注意在上面所提到最简单情形,下面的条件是满足的) 422
422 EM 算法, 常称为 GEM 算法. [注] 在应用 EM 方法时, 为了避免M-步骤中求最大值的复杂计算, 还可以采取其它的灵活的替代 方法．例如有 ECM 算法(Conditional Maximum, CM), 在多个参数的时候，例如 ( , ) J = J1 J2 的情形，如下的交 替地求条件极值的最大化的方法，也常被用来代替 M-步骤： (1) 先取 ( 1) 2 n+ J 满足 = + (( , ' | ) ( 1) ( ) 2 ( ) 1 n n n Q J J J max (( , ' ) | ) ( ) 2 ( ) 2 ' 1 n n J Q J J J ; (2) 再取 ( 1) 1 n+ J 满足 (( , ) | ( , )) ( 1) 2 ( ) 1 ( 1) 2 ( 1) 1 n+ n+ n n + Q J J J J max (( ', )| ( , )) ( 1) 2 ( ) 1 ( 1) 1 ' 1 2 + + = n n n J Q J J J J ， 称之为 ECM 算法。一般地，ECM 算法比 EM 算法达到稳定的时间长. 另有一种 ECME 算法 (混合算法)，它就是交替地使用 ECM 算法和 EM 算法． * 2. 在数据不完全时，用增补潜在数据，对参数的 Bayes 分布作估计 – Tanner–Wong 的潜变量法 ２．１ 基本想法－估计后验分布 在数据资料不全（缺失数据）, 或观测到的资料 Y 并不是状态变量时, 估计状态变量 X 的分布密度 f (J, x) 中的未知参数J , 与古典统计不同处是只能用观测到的资料Y ，这时可以通过后验分布密度 p(J | Y) 得到参数J 的 Bayes 估计．然而，后验分布 p(J | Y) 一般并不知道，就需要用观测到的资料Y 对后验分布 p(J | Y) 进行估计，这就是本段的目的． Tanner – Wong 的想法是： 在某些情形下，由观测数据Y 可以通过条件分布取样的机制，来构造某 种＂增补数据＂ (记为 Z ，称为潜变量，Latent Variable)的样本值．于是 Z 的这些样本值取自条件分布密 度 p(z | y,J) ． p(z | y) 称为预测分布．而在Y = y 已知的条件下, 潜变量的各个样本值彼此是条件 独立的． 潜变量 Z 是观测不到的, 如何选好它, 最为重要． 最简单的情形是： 观测数据Y Y Y YN , , , = 1 2 L 为 N 个时刻的历史资料，其中Yi 是一个m 维向量， 如果它的第１个分量是缺失的，而其它m -1维就 是状态变量．这时候潜变量就可取为那个缺失的一维变量． 然而，一般情形远非如此简单． 潜变量的选取是关系到能否有效地计算的关键． Tanner-Wong选取 Z 的原则是，同时满足以下两个条件（注意在上面所提到最简单情形，下面的条件是满足的）：
(1)易于模拟条件分布P(二|y9)随机数, (2)容易计算后验分布密度p(9|y,=) 下面的讨论都基于(1),(2)均能够施行 2.2未知参数的后验分布的迭代估计 旦选定潜变量以后,利用广义全概率公式,就得到 p(9|)=p9|,)(|y) p(|1)=「pe|9)p(9|y)d9 这两个关系式是设计迭代算法的基点 (迭代可能性的理论依据粗略地为:把(15.7)代入(15.6)得到 p(91y)=[ p(913,2)P(=ly, 9')d=]p(9'lyydg =JK(99m(91y)9 这说明Pp(9|y)是积分核K(9,9)的不变函数.积分核K(9,9)依赖于观测值y,但是观测值y是 固定的,所以没有将它计入记号.只要对积分核K(9,9)作适当的假定,就可以使迭代算法收敛到积分 核的不变函数) 由于积分核是不知道的,我们就不可能直接利用以下的迭代算法 gn(9)=「K(9.9gn(9d9 近似其不变函数, Tanner-Wong提出用 Monte carlo迭代方法给出p(9|y)的一个估计 Tanner-Wong 设计的迭代算法是,利用按预测分布的多个独立取样,得到潜变量的多个样本值,通过它们由后验分布 p(9|,=)们更新未知参数9的后验密度p(9|j)迭代估计.注意(5.6说明p(9|是 p(9|y,二)关于预测分布密度p(=|y)的数学期望。因此如果得到了潜变量的多个样本值 (1)→(m) 就可以作P(9|y)的如下的估计 p(9|y)=-∑p(9|y=) (15.8) 对此还需要计算潜变量加入条件的后验分布p(9|y,=) 于是,估计未知参数9的后验密度p(9|y)就归结为以下的迭代程序 (1)置初值P0(9|y)
423 （１）易于模拟条件分布 p(z | y,J) 随机数， （２）容易计算后验分布密度 p(J | y,z) ． 下面的讨论都基于（１），（２）均能够施行． 2．2 未知参数的后验分布的迭代估计 一旦选定潜变量以后, 利用广义全概率公式, 就得到 p(J | y) p(J | y,z) p(z | y)dz ò = (15. 6) 和 * * * p(z | y) p(z | y,J ) p(J | y)dJ ò = ． (15. 7) 这两个关系式是设计迭代算法的基点． （迭代可能性的理论依据粗略地为： 把 (15. 7)代入(15. 6)得到 * * * p(J | y) [ p(J | y,z) p(z | y,J )d z]p(J | y)dJ ò ò = * * * K(J,J )p(J | y )dJ ò = 记为 ． 这说明 p(J | y) 是积分核 ( , ) * K J J 的不变函数．积分核 ( , ) * K J J 依赖于观测值 y ，但是观测值 y 是 固定的，所以没有将它计入记号. 只要对积分核 ( , ) * K J J 作适当的假定， 就可以使迭代算法收敛到积分 核的不变函数）． 由于积分核是不知道的，我们就不可能直接利用以下的迭代算法 * * n * n g J K(J,J )g (J )dJ ò +1 ( ) = 近似其不变函数． Tanner – Wong 提出用Monte Carlo 迭代方法给出 p(J | y) 的一个估计．Tanner – Wong 设计的迭代算法是，利用按预测分布的多个独立取样，得到潜变量的多个样本值， 通过它们由后验分布 ( | , ) (i) p J y z 们更新未知参数J 的后验密度 p(J | y) 迭代估计．注 意(15. 6)说明 p(J | y) 是 p(J | y,z) 关于预测分布密度 p(z | y) 的数学期望。 因 此如果得到了潜变量的多个样本值 (1) ( ) , , m z L z , 就可以作 p(J | y) 的如下的估计 ( | , ) 1 ( | ) ( ) 1 ^ i m i p y z m p J y å J = = ． (15. 8) 对此还需要计算潜变量加入条件的后验分布 p(J | y,z) ． 于是, 估计未知参数J 的后验密度 p(J | y) 就归结为以下的迭代程序: （1） 置初值 ( | ) 0 p J y ;