Table Of ContentFORWARD-REVERSE EM ALGORITHM FOR MARKOV CHAINS:
CONVERGENCE AND NUMERICAL ANALYSIS
CHRISTIANBAYER,HILMARMAI,ANDJOHNSCHOENMAKERS
5
1
0 Abstract. Wedevelopaforward-reverseEM(FREM)algorithmforestimatingparame-
tersthatdeterminethedynamicsofadiscretetimeMarkovchainevolvingthroughacertain
2
measurablestatespace. AsakeytoolfortheconstructionoftheFREMmethodwede-
n velopforward-reverserepresentationsforMarkovchainsconditionedonacertainterminal
a state. TheserepresentationsmaybeconsideredasanextensionoftheearlierworkBayer
J andSchoenmakers[2013]onconditionaldiffusions.Weproofalmostsureconvergenceof
8 ouralgorithmforaMarkovchainmodelwithcurvedexponentialfamilystructure.Onthe
2 numericalsidewegiveacomplexityanalysisoftheforward-reversealgorithmbyderiving
itsexpectedcost.Twoapplicationexamplesarediscusstodemonstratethescopeofpossi-
] bleapplicationsrangingfrommodelsbasedoncontinuoustimeprocessestodiscretetime
T
Markovchainmodels.
S
.
h
t
a
m 1. Introduction
[
The EM algorithm going back to the seminal paper Laird and Rubin [1977] is a very
1 general method for iterative computation of maximum likelihood estimates in the setting
v ofincompletedata. Thealgorithmconsistsofanexpectationstep(E-step)followedbya
1 maximizationstep(M-step)whichledtothenameEMalgorithm.
9
Duetoitsgeneralapplicabilityandrelativesimplicityithasnowadaysfounditswayinto
0
a great number of applications. These include maximum likelihood estimates of hidden
7
0 MarkovmodelsinMacDonaldandZucchini[1997],non-lineartimeseriesmodelsinChan
. andLedolter[1995]andfullinformationitemfactormodelsinMengandSchilling[1996]
1
0 togivejustaverylimitedselection.
5 Despite the simplicity of the basic idea of the algorithm its implementation in more
1 complexmodelscanberatherchallenging. Theglobalmaximizationofthelikelihoodin
:
v theM-stephasrecentlybeenaddressedsuccessfully(seee.g. MengandRubin[1993]and
i Liu and Rubin [1994]). On the other hand, when the expectation of the complete likeli-
X
hoodisnotknowninclosedformonlypartialsolutionshavebeengivenyet.Oneapproach
r
developed in Wei and Tanner [1990] uses Monte Carlo approximations of the unknown
a
expectationandwasthereforenamedMonteCarloEM(MCEM)algorithm. Asanalterna-
tiveprocedurethestochasticapproximationEMalgorithmwassuggestedinLavielleand
Moulines[1999].
Date:January29,2015.
2000MathematicsSubjectClassification. 65C05,65J20.
Keywordsandphrases. EMalgorithm,Forward-reverserepresentations,Markovchainestimation,maximum
likelihoodestimation,MonteCarlosimulation.
PartiallysupportedbytheDFGResearchCenterMatheon“MathematicsforKeyTechnologies”inBerlinand
the”FinanceandSustainableGrowth”LaboratoryoftheLouisBachelierLABEXinParis.
C.B.isgratefultoPedroVillanovaforenlighteningdiscussionsregardingtheimplementationofthealgorithm
leadingtosubstantialimprovementsinTheorem6.3.
1
2 CHRISTIANBAYER,HILMARMAI,ANDJOHNSCHOENMAKERS
Inthispaperwetakeacompletelydifferentroutebyusingaforward-reversealgorithm
(cf. Bayer and Schoenmakers [2013]) to approximate the conditional expectation of the
completedatalikelihood.InthisrespectweextendtheideafromBayerandSchoenmakers
[2013] to a Markov chain setting, which is considered an interesting contribution on its
own. Indeed, Markov chains are more general in a sense, since any diffusion monitored
atdiscretetimesyieldscanonicallyaMarkovchain,butnoteverychaincanbeembedded
(straightforwardly)intosomecontinuoustimediffusiontheotherwayaround.
The central issue is the identification of a parametric Markov chain model (X , n =
n
0,1,...) based on incomplete data, i.e. realizations of the model, given on a typically
coursegridoftimepoints,letussayn ,n ,...n .Letusassumethatthechainrunsthrough
1 2 N
Rd andthatthetransitiondensities pθ (x,y),n ≥ m,ofthechainexist(with pθ (x,y) :=
n,m n,n
δ (y)),wheretheunknownparameterθhastobedetermined. Thelog-likelihoodfunction
x
basedontheincompleteobservations(X ,...,X )isthengivenby
n1 nN
(cid:88)N−1
(1.1) l(θ,x ,...,x )= lnpθ (x ,x ),
n1 nN ni,ni+1 ni ni+1
i=0
with X = x beingtheinitialstateofthechain. Thenthestandardmethodofmaximum
n0 0
likelihoodestimationwouldsuggesttoevaluate
(cid:88)N−1
(1.2) argmaxl(θ,X)=argmax lnpθ (X ,X ).
ni,ni+1 ni ni+1
θ θ i=0
Theprobleminthisapproachliesinthefactthatusuallyonlytheone-steptransitionden-
sities pθ (x,y) are explicitly known, while any multi-step density pθ (x,y) for m > n
n,n+1 n,m
canbeexpressedasanm−n−1foldintegralofone-stepdensities. Inparticularforlarger
m−n,thesemultipleintegralsarenumericallyintractablehowever.
IntheEMapproach,wethereforeconsiderthealternativeproblem
(cid:88)N−1n(cid:88)i+1−1
(1.3) argmax Elnpθj,j+1(Xj,Xj+1),
θ i=0 j=ni
intermsofthe“missingdata”Xni+1,...,Xni+1−1,i=0,...,N−1.Assuch,betweentwosuch
consecutivetimepoints,ni andni+1 say,thechainmaybeconsideredasabridgeprocess
startinginrealizationX andendingupinrealizationX (undertheunknownparameter
ni ni+1
θ though), and so each term in (1.3) may be considered as an expected functional of the
“bridged”Markovchainstartingattimen in(datapoint)X ,conditionalonreaching(data
i ni
point)Xni+1 attimeni+1.
We will therefore develop firstly an algorithm for estimating the terms in (1.3) for a
givenparameterθ.Thisalgorithmwillbeofforward-reversetypeinthespiritoftheone
in Bayer and Schoenmakers [2013] developed for diffusion bridges. It should be noted
herethatinthelastyearstheproblemofsimulatingdiffusionbridgeshasattractedmuch
attention. Withoutpretendingtobecomplete,seeforexample,BladtandSørensen[2014],
DelyonandHu[2006],MilsteinandTretyakov[2004],Stinis[2011],Stuartetal.[2004],
Schaueretal.[2013].
Havingtheforward-reversealgorithmathand,wemayconstructanapproximatesolu-
tionto(1.3)inasequentialwaybytheEMalgorithm: Onceagenericapproximationθ is
m
constructedaftermsteps,oneestimates
(cid:88)N−1n(cid:88)i+1−1
θm+1 :=argmax (cid:98)Elnpθj,j+1(Xθjm,Xθj+m1),
θ i=0 j=ni
FORWARD-REVERSEEMALGORITHMFORMARKOVCHAINS 3
whereXθm denotestheMarkovbridgeprocessunderthetransitionlawduetoparameterθm
andeachterm
(cid:98)Elnpθ (Xθm,Xθm )
j,j+1 j j+1
representsaforward-reverseapproximationof
Elnpθ (Xθm,Xθm )
j,j+1 j j+1
asa(known)functionofθ.
ConvergencepropertiesofapproximateEMalgorithmshavedrawnconsiderablerecent
attention in the literature mainly driven by it’s success in earlier intractable estimation
problems. An overview of existing convergence results for the MCEM algorithm can be
fund in Neath [2013]. Starting from a convergence result for the forward-reverse repre-
sentation for Markov chains we prove almost sure convergence of the FREM sequence
inthesettingofcurvedexponentialfamiliesbasedontechniquesdevelopedinBayerand
Schoenmakers[2013]andFortandMoulines[2003]. Essentiallytheonlyingredientfrom
theMarkovchainmodelforthisconvergencetoholdareexponentialtailsofthetransition
densitiesandtheirderivativesthatarestraightforwardtocheckinmanyexamples.
Sincecomputationalcomplexityisalwaysanissueinestimationtechniquesthatinvolve
asimulationstep,wealsoincludeacomplexityanalysisfortheforward-reversealgorithm.
WeshowthatthealgorithmachievesanexpectedcostoftheorderO(NlogN)forsample
sizeofN forward-reversetrajectories.
WementionarecentapplicationpaperbyBayer,Moraes,Tempone,andVilanovaBayer
etal.[2015],whichfocusesonpracticalandimplementationissuesoftheforward-reverse
EMalgorithminthesettingofStochasticReactionNetworks(SRNs), i.e., ofcontinuous
timeMarkovchainswithdiscrete,butgenerallyinfinitestatespace.
The structure of the paper is as follows. In Section 2 we recapitulate and adapt the
conceptofreversedMarkovchains,initiallydevelopedinMilsteinetal.[2007]usingthe
ideasinMilsteinetal.[2004]onreverseddiffusions. Ageneralstochasticrepresentation
— involving standard (unconditional) expectations only — for expected functionals of
conditional Markov chains is constructed in Section 3. This representation allows for a
forwardreverseEMalgorithmthatisintroducedandanalyzedinSection4. InSection5
we proof almost sure convergence of the forward-reverse EM algorithm in the setting of
curvedexponentialfamilies. ImplementationandComplexityoftheFREMalgorithmare
addressedinSection6. ThepaperisconcludedwithtwoapplicationexamplesinSection
7thatdemonstratethewidescopeofourmethod.
2. RecapofforwardandreverserepresentationsforMarkovchains
Consideradiscrete-timeMarkovprocess(X ,F ), n=0,1,2,...,onaprobabilityspace
n n
(Ω,F,P)withphasespace(S,S),henceforthcalledMarkovchain. Ingeneralweassume
thatS islocallycompactandthatSistheBorelσ-algebraonS.Forexample,S =Rd ora
propersubsetofRd.LetP , n≥0,denotetheone-steptransitionprobabilitiesdefinedby
n
(2.1) Pn(x,B):=P(Xn+1 ∈ B|Xn = x), n=0,1,2,..., x∈S, B∈S.
InthecaseofanautonomousMarkovchainalltheone-steptransitionprobabilitiescoincide
andareequaltoP:= P = P =···.
0 1
Let Xn,x, m ≥ n, be atrajectory of theMarkov chain whichis at step n in thepoint x,
m
i.e.,Xn,x = x.Themulti-steptransitionprobabilitiesP arethendefinedby
n n,m
P (x,B):=P(Xn,x ∈ B), x∈S, B∈S, m≥n.
n,m m
4 CHRISTIANBAYER,HILMARMAI,ANDJOHNSCHOENMAKERS
Duetothesedefinitions,Pn,n(x,B) = δx(B) = 1B(x)(Diracmeasure),Pn = Pn,n+1,andthe
Chapman-Kolmogorovequationhasthefollowingform:
(cid:90)
(2.2) P (x,B)= P (x,dy)P (y,B), x∈S, B∈S, n≤k≤m.
n,m n,k k,m
LetusfixN >0andconsiderfor0≤n≤ N thefunction
(cid:90)
(2.3) u (x):= P (x,dy)f(y)=E f(Xn,x),
n n,N N
where f is S-measurable and such that the mathematical expectation in (2.3) exists; for
example, f isbounded. BytheMarkovpropertywehavefor0≤n< N :
u (x)=E f(Xn,x)=E f(Xn+1,Xnn+,x1)
n N N
=EEFn+1f(Xn+1,Xnn+,x1)=EEXnn+,x1f(Xn+1,Xnn+,x1)
N N
(cid:90)
=Eun+1(Xnn+,x1)= un+1(y)Pn(x,dy).
Thus,u (x)satisfiesthefollowingdiscreteintegralCauchyproblem
n
(cid:90)
(2.4) un(x)= un+1(y)Pn(x,dy), n< N,
(2.5) u (x)= f(x),
N
and(2.3)isaforwardprobabilisticrepresentationofitssolution. Infact,theprobabilistic
representation(2.3)canbeusedforsimulatingthesolutionof(2.4)-(2.5)byMonteCarlo.
For our purpose, reverse probabilistic representations we need a somewhat more general
versionoftheaboveresult.
Theorem 2.1 (cf. Milstein et al. [2007]). Let P be the one-step transition density of a
n
Markovchain X asin(2.1)andletthefunction f :S → Rbemeasurableandbounded.
Letfurtherϕ :S×S →Rbeameasurableandboundedfunctionsforn=0,1,2,...Then,
n
thesolutionoftheproblem
(cid:90)
(2.6) wn(x)= wn+1(z)ϕn(x,z)Pn(x,dz), n< N,
(2.7) w (x)= f(x)
N
hasthefollowingprobabilisticrepresentation:
(cid:104) (cid:105)
(2.8) w (x)=E f(Xn,x)Xn,x,1 ,
n N N
where(X,X)isanextendedMarkovchaininwhichXisgovernedbytheequations
Xn,x,γ =Xn,x,γϕ (Xn,x,Xn,x), Xn,x,γ =γ,
k+1 k k k k+1 n
wheren≤k< N.
FORWARD-REVERSEEMALGORITHMFORMARKOVCHAINS 5
Proof. NotethatXn,x,γ =γXn,x,1.Thus,forn< N,(2.8)maybewrittenas
k k
(cid:20) (cid:21)
w (x)=E f(Xn+1,Xnn+,x1)Xn+1,Xnn+,x1,Xnn,+x1,1
n N N
(cid:20) (cid:21)
=EXnn,+x1,1E(Xnn+,x1,Xnn,+x1,1) f(XNn+1,Xnn+,x1)XnN+1,Xnn+,x1,1
(cid:104) (cid:105)
=E Xnn,+x1,1wn+1(Xnn+,x1)
(cid:104) (cid:105)
=E ϕn(x,Xnn+,x1)wn+1(Xnn+,x1)
(cid:90)
= wn+1(z)ϕn(x,z)Pn(x,dz),
and(2.7)istriviallyfulfilledforn= N. (cid:3)
(cid:16) (cid:17)
2.1. Reverse probabilistic representations. We henceforth take (S,S) = Rd,B(Rd)
andassumethatthetransitionprobabilitiesP (x,dy)havedensitiesp (x,y)withrespect
n,m n,m
to the Lebesgue measure on (S,S). We note however that without any problem one may
considermoregeneralstatespacesequippedwithsomereferencemeasure,andtransition
probabilitiesabsolutelycontinuoustowithrespecttoit. Therepresentation(2.3)canthus
bewrittenintheform
(cid:90)
(2.9) I(f):=E f(Xn,x)= p (x,y)f(y)dy, 0≤n≤ N.
N n,N
Lettheinitialvalueξ ofthechain X atmomentnberandomwithdensityg(x).Consider
thefunctional
(cid:90) (cid:90)
(2.10) I(g, f)= g(x)p (x,y)f(y)dxdy=Ef(Xn,ξ).
n,N N
Formally, by taking for g a δ-function we obtain (2.9) again, and by taking f to be a δ-
functionweobtaintheintegral
(cid:90)
(2.11) J(g):= g(x)p (x,y)dx.
n,N
We now propose suitable (reverse) probabilistic representations for J(g), where g is an
arbitrarytestfunction(notnecessarilyadensity). Forthiswearegoingtoconstructaclass
of reverse Markov chains that allow for a probabilistic representation for the solution of
(2.11).
Letusfixanumber N ∈ Nandconsiderfor0 ≤ m < N,functionsψm :S ×S → R+
suchthatforeachmandythefunction
p (·,y)
(2.12) q (y,·):= N−m−1 , 0≤m< N,
m ψ (y,·)
m
is a density on S. For example, one could take ψ independent of the second argument,
m
andthenobviously
(cid:90)
(2.13) ψ (y)= p (z,y)dz.
m N−m−1
Wenowintroducea“reverse”processes(Yy,Yy) bythesystem
m m 0≤m≤N
(cid:12)
P(Ymy+1 ∈dz(cid:48)(cid:12)(cid:12) Ymy =z)=qm(z,z(cid:48))dz(cid:48),
(2.14) Yy =Yyψ (Yy,Yy ),
m+1 m m m m+1
Yy :=Y0,y :=y, Yy :=Y0,y,1 :=1, 0≤m< N,
0 0 0 0
6 CHRISTIANBAYER,HILMARMAI,ANDJOHNSCHOENMAKERS
henceYy isgovernedbytheone-steptransitionprobabilitiesQ (z,dz(cid:48)):=q (z,z(cid:48))dz(cid:48) (i.e.
m m
Q insteadofP ).
m m
Theorem2.2. Foranyn,0≤n≤ N,(2.11)hasthefollowingprobabilisticrepresentation.
(cid:90)
(cid:104) (cid:105)
g(x)p (x,y)dx=E g(Yy )Yy ,
n,N N−n N−n
where g is an arbitrary test function (a ”density” p has to be interpreted as a Dirac
m,m
distributionorδ-function).
Proof. From the Chapman - Kolmogorov equation (2.2) we obtain straightforwardly the
Chapman-Kolmogorovequationfordensities,
(cid:90)
(2.15) p (x,y)= p (x,z)p (z,y)dz, x,y∈S, n≤k≤m.
n,m n,k k,m
Letusnowfixn,n< N(forn= Nthestatementistrivial)also,andintroducethefunctions
(cid:90)
(2.16) v (y):= g(x)p (x,y)dx, n≤k≤ N.
k n,k
From(2.15)weget
(cid:90)
(2.17) v (y)= v (z)p (z,y)dz, n<k≤ N,
k k−1 k−1
v (y)=g(y),
n
where p := p denote the one-step densities. For n < k ≤ N we now consider a
k−1 k−1,k
“reversed”timevariablem= N+n−kandwritewith(cid:101)vm(y):=vN+n−m(y)and(2.12)system
(2.17)intheform
(cid:90)
(2.18) (cid:101)vm(y)= (cid:101)vm+1(z)ψm−n(y,z)qm−n(y,z)dz, n≤m< N,
(cid:101)vN(y)=g(y).
Letuswrite(2.18)inaslightlydifferentform,
(cid:90)
(cid:101)vm(y)= (cid:101)vm+1(z)ψ(mn)(y,z)q(mn)(y,z)dz, n≤m< N,
(cid:101)vN(y)=g(y)
withψ(n) := ψ andq(n) := q .ViaTheorem2.1wenextobtainaprobabilisticrepre-
m m−n m m−n
sentationoftheform(2.8)forthesolutionofproblem(2.18),hence(2.11)orJ(g).Indeed,
(cid:16) (cid:17)
bytakinginTheorem2.1insteadof X aMarkovchain Y(n),y ,whereY(n),y isgov-
m n≤m≤N
ernedbytheone-steptransitionprobabilitiesQ(n)(z,dz(cid:48)) :=q(n)(z,z(cid:48))dz(cid:48),n ≤ m < N,with
(cid:16) m (cid:17) m
initialconditionY(n),y =y,andconstructing Y(n),y accordingto
n m n≤m≤N
(2.19) Y(n),y =Y(n),yψ(n)(Y(n),y,Y(n),y), Y(n),y =1, n≤m< N,
m+1 m m m m+1 n
itfollowsbyTheorem 2.1that
(cid:104) (cid:105)
(2.20) J(g)=(cid:101)vn(y)=vN(y)=E g(YN(n),y)Y(Nn),y .
Itremainstoseethat
(cid:104) (cid:105) (cid:104) (cid:105)
E g(Y(n),y)Y(n),y =E g(Yy )Yy
N N N−n N−n
FORWARD-REVERSEEMALGORITHMFORMARKOVCHAINS 7
which follows from the fact that initial values and the one step transition probabilities of
theprocesses
(cid:16) (cid:17) (cid:16) (cid:17)
Y(n),y,Y(n),y and Yy,Yy
n+i n+i i=0,...,N−n i i i=0,...,N−n
coincide. (cid:3)
It should be stressed that, in contrast to a corresponding theorem in Milstein et al.
[2007], Theorem 2.2 provides a family of probabilistic representations indexed by n =
1,...,N,thatinvolvesonlyonecommonreverseprocessYy.InTheorem2.2N wasfixed
but,whendifferentN areinplay,wewilldenotethembyYy;N.Itturnsoutthatthisexten-
sionoftherelatedresultinMilsteinetal.[2007]iscrucialforderivingprobabilisticrepre-
sentationsforconditionalMarkovchainsbelow(cf. BayerandSchoenmakers[2013]).
3. Simulationofconditionalexpectationsviaforward-reverserepresentations
InthissectionwedescribeforaMarkovChain(2.1)anefficientprocedureforestimat-
ingthefinaldistributionsofachainX = (Xn)n=0,...,N conditioned,orpinned,onaterminal
stateX .Morespecifically,forsomegiven(unconditional)processXweaimatsimulation
N
ofthefunctional
(cid:104) (cid:12) (cid:105)
(3.1) E g(Xm1,...,Xmr)(cid:12)(cid:12) XN =y, X0 = x ,
where 0 ≤ m < m < ··· < m < N (hence r < N), g is an arbitrarily given suitable
1 2 r
test function, and x,y ∈ Rd are given states. The procedure proposed below is in fact an
extension of the method developed in Bayer and Schoenmakers [2013] to discrete time
Markov chains. We note that similar techniques as in Bayer and Schoenmakers [2013]
alsoallowustotreatthemoregeneralproblem
(cid:104) (cid:12) (cid:105)
E g(Xm1,...,Xmr)(cid:12)(cid:12) XN ∈ A, X0 = x
forsuitablesetsA⊂Rd.
3.1. Forward-reverserepresentationsofconditionalexpectations. Letusconsiderthe
problem (3.1) for fixed x,y ∈ Rd (i.e. A = {y}). We firstly state the following central
theorem.
Theorem3.1. GivenagridD :={0≤n∗ <n <···<n =: N},itholdsthat
l 1 l
(cid:104) (cid:105)
E f(Yy;nl ,Yy;nl ,...,Yy;nl )Yy;nl
nl−n0 nl−n1 nl−nl−1 nl−n0
(cid:90) (cid:89)l
= f(y ,y ,...,y ) p (y ,y)dy
0 1 l−1 ni−1,ni i−1 i i−1
Rd×L i=1
withy :=yandn :=n∗.
l 0
Proof. Withoutlossofgenerality,weassumeinthisproofthatthegridsatisfiesn −n =
i i−1
1,i=1,...,l. Indeed,extend f :Rd×l →Rtoafunction f˜:Rd×(N−n∗) →Rsuchthat
(cid:16) (cid:17) (cid:16) (cid:17)
f˜ Yy;N ,Yy;N ,...,Yy;N,Yy;N = f Yy;N ,Yy;N ,...,Yy;N .
N−n∗ N−n∗−1 2 1 nl−n0 nl−n1 nl−nl−1
Then,re-expressingthetransitiondensitiesp intermsoftheone-steptransitiondensi-
ni−1,ni
ties p usingChapman-Kolmogorov,weseethatthestatementofthetheoremisequivalent
i
8 CHRISTIANBAYER,HILMARMAI,ANDJOHNSCHOENMAKERS
to
(cid:104) (cid:16) (cid:17) (cid:105)
(3.2) E f˜ Yy;N ,Yy;N ,...,Yy;N Yy;N
N−n∗ N−n∗−1 1 N−n∗
(cid:90) (cid:89)N
= f˜(y ,...,y ) p (y ,y)dy
n∗ N−1 i−1 i−1 i i−1
Rd×(N−n∗) i=n∗+1
withy ≡y. Infact,weshallprovethat
N
(cid:104) (cid:16) (cid:17) (cid:105) (cid:90) (cid:89)N
(3.3) E f Yy;N,...,Yy;N Yy;N = f (y ,...,y ) p (y ,y)dy
p p 1 p p N−p N−1 i−1 i−1 i i−1
i=N−p+1
forany1 ≤ p ≤ N−n∗ forany(e.g.,boundedmeasurable)function f : Rd×p → R. (3.3)
p
givestheformulafromthestatementofthetheoremfor p = N −n∗ with f beingthe
N−n∗
function f˜fromabove. Weprove(3.3)byinductionon p. For p = 1, thisboilsdownto
Theorem2.2withn= N−1.
Forthestepfrom p−1to p,wenotethatbydefinition
(cid:16) (cid:17)
Yy;N =Yy;Nψ Yy;N,Yy;N ,
p p−1 p−1 p−1 p
withψ (y,·)q (y,·)= p (·,y)= p (·,y)by(2.12). Hence,wehave
p−1 p−1 N−(p−1)−1 N−p
(cid:104) (cid:105) (cid:104) (cid:105)
E f (Yy;N,Yy;N,...,Yy;N)Yy;N =E f (Yy;N,Yy;N,...,Yy;N)Yy;Nψ (Yy;N,Yy;N)
p p p−1 1 p p p p−1 1 p−1 p−1 p−1 p
(cid:104) (cid:16) (cid:17) (cid:105)
=E g Yy;N,...,Yy;N Yy;N ,
p−1 1 p−1
with
g(z ,...,z )≡E(cid:20)f (cid:16)Yy;N,Yy;N,...,Yy;N(cid:17)ψ (Yy;N,Yy;N)(cid:12)(cid:12)(cid:12)Yy;N =z ,...,Yy;N =z (cid:21)
p−1 1 p p p−1 1 p−1 p−1 p (cid:12) p−1 p−1 1 1
(cid:90)
= f (z,z ,...,z )p (z,z )dz.
p p−1 1 N−p p−1
Applyingtheinductionhypothesisfor f =g,weobtain
p−1
(cid:104) (cid:105) (cid:104) (cid:16) (cid:17) (cid:105)
E f (Yy;N,Yy;N,...,Yy;N)Yy;N =E g Yy;N,...,Yy;N Yy;N
p p p−1 1 p p−1 1 p−1
(cid:90) (cid:89)N
= g(yN−p+1,...,yN−1) pi−1(yi−1,yi)dyi−1
i=N−p+2
(cid:90) (cid:89)N
= f (y ,...,y ) p (y ,y)dy . (cid:3)
p N−p N−1 i−1 i−1 i i−1
i=N−p+1
Wenowconsideranextendedintegersequence
0<m <···<m =n∗ =n <n <···<n = N,
1 k 0 1 l
andakernelK oftheform
(cid:15)
K (u):=(cid:15)−dK(u/(cid:15)), y∈Rd,
(cid:15)
(cid:82)
with K being integrable on Rd and K(u)du = 1. Formally K converges to the delta
Rd (cid:15)
function δ on Rd (in distribution sense) as (cid:15) ↓ 0. We then have the following stochastic
0
representation(involvingstandardexpectationsonly)for(3.1) withni =mk+i,i=0,...,l=
r−k+1.
FORWARD-REVERSEEMALGORITHMFORMARKOVCHAINS 9
(cid:16) (cid:17)
Theorem 3.2. Let the chain (Y,Y) := Yy;N,Yy;N be given by (2.14), and the modified
(cid:0) (cid:1)
integersequence (cid:98)n· bedefinedby
(3.4) (cid:98)ni :=nl−nl−i, i=1,...,l.
Itthenholds
(cid:104) (cid:12) (cid:105)
E g(Xm1,...,Xmr)(cid:12)(cid:12) Xm0 = x, XN =y
(cid:104) (cid:12) (cid:105)
=E g(Xm1,...,Xmk−1,Xn0∗,x,Xn1,...,Xnl−1)(cid:12)(cid:12) Xm0 = x, XN =y
(cid:104) (cid:16) (cid:17) (cid:16) (cid:17) (cid:105)
E g Xm0,x,...,Xm0,x,Xm0,x,Yy;N,...,Yy;N K Yy;N −Xm0,x Yy;N
(3.5) =lim m1 mk−1 (cid:104)n∗ (cid:16) (cid:98)nl−1 (cid:98)n(cid:17)1 (cid:15)(cid:105) (cid:98)nl n∗ (cid:98)nl .
(cid:15)→0 E K Yy;N −Xm0,x Yy;N
(cid:15) (cid:98)nl n∗ (cid:98)nl
Proof. TheproofisanalogoustothecorrespondingoneinBayerandSchoenmakers[2013].
Asaroughsketch,applyTheorem3.1to
f(X0,x,...,X0,x,y ,y ,...,y ):=g(X0,x,...,X0,x,y ,...,y )K (y −X0,x),
m1 n∗ 0 1 l−1 m1 n∗ 1 l−1 (cid:15) 0 n∗
conditionalonX0,x,...,X0,x,send(cid:15) →0,anddividetheresultby
m1 n∗
p (x,y)=limE(cid:2)K (cid:0)Y −X (cid:1)Y (cid:3). (cid:3)
0,N (cid:15)→0 (cid:15) (cid:98)nl n∗ (cid:98)nl
Forward-Reversealgorithm. GivenTheorem3.2thecorrespondingforward-reverseMonte
Carloestimatorfor(3.5)suggestsitself: Samplei.i.d. copiesX0,x,(1),...,X0,x,(M) ofthepro-
(cid:16) (cid:17) (cid:16) (cid:17)
cess X0,x and, independently, i.i.d. copies Yy;N,(1),Yy;N,(M(cid:101)) ,..., Yy;N,(1),Yy;N,(M(cid:101)) of the
(cid:16) (cid:17)
process Yy;N,Yy;N . Take for K a second order kernel, take for simplicity M = M(cid:101), and
choose a bandwidth (cid:15) ∼ M−1/d if d ≤ 4, or (cid:15) ∼ M−2/(4+d) if d ≥ 4. By next replac-
M M
ingtheexpectationsinthenumeratoranddenominatorof(3.5)bytheirrespectiveMonte
Carlo estimates involving double sums, one ends up with an estimator with Root-Mean-
SquareerrorO(M−1/2)inthecased ≤4andO(M−4/(4+d))inthecased >4(cf. Bayerand
Schoenmakers[2013]fordetails).
4. Theforward-reverseEMalgorithm
Let us now formulate the forward-reverse EM (FREM) algorithm in the setting of the
missingdataproblem. Supposethattheparameterθ ∈ Θ ⊂ Rs andthattheMarkovchain
X = (X ,n ∈ N)hasstatespaceRd. Assumingthatthetransitiondensities p of X exist
n n,k
forn,k∈Nthefull-datalog-likelihoodthenreads
n(cid:88)N−1
(4.1) lc(θ,x)= logpθi,i+1(xi,xi+1), x∈RnN+1.
i=0
In the missing data problem only partial observations X ,...,X are available for 0 =
n0 nN
n < n < ... < n with log-likelihood function l given in (1.1) that is intractable in
0 1 N
mostcases. Instead, themaximizationoflinθ hastobereplacedbyatwostepiterative
procedure,theEMalgorithm.
E-step: In the m-th step evaluate the conditional expectation of the complete data
log-likelihood
Q(θ,θ ,x):=E [l (θ,X ,...,X )|X = x ,...,X = x ], x∈RN+1.
m θm c 0 nN n0 n0 nN nN
M-step: Updatetheparameterby
θm+1 =argmaxQ(θ,θm,X).
θ
10 CHRISTIANBAYER,HILMARMAI,ANDJOHNSCHOENMAKERS
Since in many Markov chain models the E-step is intractable in this form, we propose a
forward-reverseapproximationfor theexpectationofthetransition densitiesevaluatedat
theobservations.
FRE-step: Evaluate
(4.2) Q (θ,θ ,X):=EFR[l (θ,X)|X ,i=0,...,N],
m m θm c ni
whereEFRdenotesaforward-reverseapproximationoftheconditionalexpectation
θm
undertheparameterθ .
m
After this FR E-step is computed the M-step remains unchanged. This FREM algorithm
gives a random sequence (θ ) that under certain conditions given in the next section
n n≥0
converges to stationary points of the likelihood function. To assure a.s. boundedness of
thissequenceweapplyastabilizationtechniqueintroduceinChenetal.[1988].
ThestableFREMalgorithm. Let K ⊂ Θform ∈ Nbeasequenceofcompactsetssuch
m
that
(cid:91)
(4.3) Km (cid:40) Km+1 and Θ= Km
m∈N
for all m ∈ N. We define the stable FREM algorithm by checking if θ after the m-th
m
maximizationstepliesinK andresetingthealgorithmotherwise. Chooseastartingvalue
m
θ ∈ K andlet p forn∈N, p :=0,countthenumberofresets.
0 0 n 0
stableM-step:
(4.4) θm+1 =argmaxQm(θ,θm,X)and pn+1 = pn, if argmaxQm(θ,θm,X)∈ Km,
θ θ
(4.5) θm+1 =θ0and pn+1 = pn+1, if argmaxQm(θ,θm,X)(cid:60) Km.
θ
Wewillshowinthenextsectionthatunderweakassumptionthenumberofresets p stays
n
a.s. finite. OurstableFREMalgorithmconsistsnowofiterativelyrepeatingtheFRE-step
andthestableM-step.
5. AlmostsureconvergenceoftheFREMalgorithm
InthissectionweprovealmostsureconvergenceofthestableFREMalgorithmunder
theassumptionthatthecompletedatalikelihoodisfromacurvedexponentialfamily. Our
proofismainlybasedonresultsfromBayerandSchoenmakers[2013],FortandMoulines
[2003] and the classical framework for the EM algorithm introduced in Laird and Rubin
[1977]andLange[1995].
5.1. Model setting. Suppose that φ : Θ → R, ψ : Θ → Rq and S : R(N+1)d → Rq are
continuousfunctions. Wemakethestructuralassumptionthatthefulldatalog-likelihood
isoftheform
(5.1) l(θ,x ,...,x )=φ(θ)+(cid:104)S(x ,...,x ),ψ(θ)(cid:105),
0 N 0 N
i.e. l is from a curved exponential family. In order to proof convergence we need the
followingpropertiestobefulfilledthatnaturallyholdinmanypopularmodels. InSection
6wegiveseveralpracticalexamplesthatfallintothissetting.
Assumption5.1. (1) Thereexistsacontinuousfunctionθ¯ :Rq →Θsuchthatl(θ¯(s),s)=
sup l(θ,s)foralls∈R(N+1)d.
θ∈Θ
(2) TheincompletedatalikelihoodLiscontinuousinθ,andthelevelsets{θ∈Θ|L(θ,x)≥
C}arecompactforanyC >0andallx.
