Reviving the Two-state Markov Chain Approach (Technical Report)

0.84 MB·English

by  Andrzej Mizera

Reviving the Two-state Markov Chain Approach
(Technical Report)
AndrzejMizera1,JunPang1,2,andQixiaYuan1(cid:63)
1 UniversityofLuxembourg,FSTC,Luxembourg
2 UniversityofLuxembourg,SnT,Luxembourg
[email protected]
Abstract. ProbabilisticBooleannetworks(PBNs)isawell-establishedcompu-
6
tationalframeworkformodellingbiologicalsystems.Thesteady-statedynamics
1
of PBNs is of crucial importance in the study of such systems. However, for
0
largePBNs,whichoftenariseinsystemsbiology,obtainingthesteady-statedis-
2
tributionposesasignificantchallenge.Infact,statisticalmethodsforsteady-state

ct approximationaretheonlyviablemeanswhendealingwithlargenetworks.In
O thispaper,werevivethetwo-stateMarkovchainapproachpresentedintheliter-
ature.Wefirstidentifyaproblemofgeneratingbiasedresults,duetothesizeof

5 theinitialsamplewithwhichtheapproachneedstostartandweproposeafew
2 heuristicstoavoidsuchapitfall.Second,weconductanextensiveexperimental

comparisonofthetwo-stateMarkovchainapproachandanotherapproachbased
]
E ontheSkartmethodandweshowthatstatisticallythetwo-stateMarkovchain
hasabetterperformance.Finally,weapplythisapproachtoalargePBNmodel
C
ofapoptosisinhepatocytes.
.
s
c
[
1 Introduction

4
v
Systems biology aims to study biological systems from a holistic perspective, with
9
thegoaltoprovideacomprehensive,system-levelunderstandingofcellularbehaviour.
7
7 Proper functioning of a living cell requires a finely-tuned and orchestrated interplay
1 of many complex processes. Complex interactions within a biological system lead to
0 emergentpropertieswhicharecrucialforsustaininglife.Therefore,understandingthe
.
1 machineryofliferequirestheuseofholisticapproacheswhichenablethestudyofasys-
0 temasawhole,incontrasttothereductionistapproach.Computationalmodellingplays
5
a prominent role in the field of systems biology. Construction and analysis of a com-
1
putational model for some biological process enables the systematisation of available
:
v biologicalknowledge,identificationofmissingbiologicalinformation,providesformal
i
X meansforunderstandingandreasoningabouttheconcertedinterplaybetweendifferent
partsofthemodel,finallyrevealsdirectionsforfutureexperimentalworkwhichcould
r
a providedataforbetterunderstandingoftheprocessunderstudy.
Unfortunately, computational modelling of biological processes that take place in
alivingcellposessignificantchallengeswithrespecttothesizeofthestate-spacethat
(cid:63)Supported by the National Research Fund, Luxembourg (grant 7814267) and partially by
GoogleSummerofCode2014.
2 Mizera,PangandYuan
needstobeconsidered.Modellingofcertainpartsofcellularmachinerysuchasgene
regulatory networks (GRNs) or signal transduction pathways often leads to dynami-
cal models characterised by huge state-spaces of sizes that surpass the sizes of any
human-designed systems by orders of magnitude. Therefore, profound understanding
ofbiologicalprocessesasksforthedevelopmentofnewmethodsandapproachesthat
wouldprovidemeansforformalanalysisandreasoningaboutsuchhugesystems.
Inthisstudyweconcentrateontheanalysisofthesteady-statedynamicsofbiologi-
calprocesses,inparticularGRNs,modelledasdiscrete-timeMarkovchains(DTMCs).
This is the case, for example, when the biological system under study is cast into the
mathematical/computationalframeworkofprobabilisticBooleannetworks(PBNs)[1,2].
Intheseorotherdiscrete-timemodels,e.g.,dynamicBayesiannetworks,thereal(con-
sidered as continuous) time is not modelled. Instead, the evolution of the system is
abstractedasasequenceofconsecutiveevents.Thesecoarse-grainedmodelshavebeen
successfullyappliedinmanysystemsbiologystudiesandprovedtheirpredictivepower[3].
In fact, for the study of large regulatory systems they remain the only reasonable so-
lution. Extrapolating the ordinary differential equations model of a single elementary
building block of the network (e.g., a gene) to the whole large system would result
in a prohibitively complex model. However, moving towards a higher-level descrip-
tionbyignoringthemoleculardetailsallowstograspthesystem-levelbehaviourofthe
network [4]. In consequence, these coarse-grained formalisms are broadly applied in
systems biology studies of systems where the predictions of exact reaction times are
notofmaininterest.Forexample,thisisthecaseinthestudyofdynamicalattractors
ofaregulatorynetwork,whichseemtodependonthecircuitwiringratherthankinetic
constants (such as production rates or interaction rates) [5]. In this sense modelling
biologicalsystemswithmoreabstract,high-levelviewformalismshascertainunques-
tionableadvantages.
Oneofthekeyaspectsintheanalysisofsuchdynamicsystemsisthecomprehen-
sionoftheirsteady-state(long-run)behaviour.Forexample,attractorsofsuchsystems
werehypothesisedtocharacterisecellularphenotypes[6].Anothercomplementarycon-
jecture is that attractors correspond to functional cellular states such as proliferation,
apoptosis, or differentiation [7]. These interpretations may cast new light on the un-
derstandingofcellularhomeostasisandcancerprogression[1].Inthiswork,wefocus
onthecomputationofsteady-stateprobabilitieswhicharecrucialforthedetermination
of long-run influences and sensitivities. These are measures that quantify the impact
of genes on other genes, considered collectively or individually, and that enable the
identificationofelementswithhighestimpact.Inthiswaytheyprovideinsightintothe
controlmechanismsofthenetwork.
Sofarthehuge-statespace,whichoftencharacterisesdynamicalmodelsofbiolog-
icalsystems,temperstheapplicationoftheabovementionedtechniquesintheanalysis
of realistic biological systems. In fact, approximations with the use of Markov chain
MonteCarlo(MCMC)techniquesaretheonlyviablesolutiontothisproblem[8].How-
ever, due to the difficulties with the assessment of the convergence rate to the steady-
statedistribution(see,e.g.,[9]),certaincareisrequiredwhenapplyingthesemethods
in practice. A number of statistical methods exists, which allow to empirically deter-
minewhentostopthesimulationandoutputestimates.Weemployinourstudy:(1)the
RevivingtheTwo-stateMarkovChainApproach 3
two-state Markov chain approach [10] and (2) the Skart batch-means method of [11].
The two-state Markov chain approach was introduced in 1992 by Raftery and Lewis;
and Shmulevich et al. [8] proposed its application to the analysis of PBNs in 2003.
However, to the best of our knowledge, since then it has not been widely applied for
theanalysisoflargePBNs.Inthispaper,weaimtoreviveitsusageforapproximating
steady-state probabilities of large PBNs, which often arise in computational systems
biologyasmodelsoffull-sizegeneticnetworks.Weidentifyaproblemconcernedwith
the choice of the initial sample size in the two-state Markov chain approach. As we
show, an unconscious choice may lead to biased results. We propose a few heuristics
to avoid this problem. By extensive experiments, we show that the two-state Markov
chain approach outperforms the Skart method in most cases, where the batch-means
Skartmethodisconsideredthecurrentstate-of-the-artapproach.Inthiswayweshow
that the two-state Markov chain is often the optimal choice for the analysis of large
PBNmodelsofbiologicalsystems.
Structureofthepaper.AfterpresentingsomepreliminariesinSection2,wedescribe
thetwo-stateMarkovchainapproachinSection3andidentifyaproblemofgenerating
biased results, due to the size of the initial sample with which the approach needs to
start. We then propose some heuristics for the approach to avoid unfortunate initiali-
sations. We perform an extensive evaluation and comparison of the two-state Markov
chainapproachandtheSkartmethodinSection4onalargenumberofrandomlygener-
atedPBNs.Inmostcases,thetwo-stateMarkovchainapproachseemstohaveabetter
performance in terms of computational cost. Finally, we apply the two-state Markov
chain approach to study a large PBN model of apoptosis in hepatocytes consisting of
91nodesinSection5.Wecomputethesteady-stateinfluencesandlong-runsensitivi-
tiesandconfirmpreviouslyformulatedhypothesis.Weconcludeourpaperwithsome
discussionsinSection6.
2 Preliminaries
2.1 Finitediscrete-timeMarkovchains(DTMCs)
LetSbeafinitesetofstates.A(first-order)discrete-timeMarkovchainisanS-valued
stochasticprocess Xt t Nwiththepropertythatthenextstateisindependentofthepast
states given the pr{ese}nt∈state. Formally, P(X =s X =s,X =s ,...,X =
t+1 t+1 t t t 1 t 1 0
s )=P(X =s X =s)foralls ,s,...,s S.He|re,wecons−idertim−e-homogenous
0 t+1 t+1 t t t+1 t 0
Markovchains,i.e.,|chainswhereP(X =s X∈=s),denotedP ,isindependentof
t+1 (cid:48) t s,s
t for any states s,s S. The transition matri|x P=(P ) sat(cid:48)isfies P (cid:62)0 and
(cid:48) s,s s,s S s,s
∑s SPs,s =1forall∈s S.Wedenotebyπ aprobability(cid:48)dist(cid:48)r∈ibutiononS.If(cid:48)π =πP,
the(cid:48)∈nπ is(cid:48)astationaryd∈istributionoftheDTMC(alsoreferredtoastheinvariantdistri-
bution).Apathoflengthnisasequences s s suchthatP >0and
1→ 2→···→ n si,si+1
s Sfori 1,2,...,n .Stateq Sisreachablefromstate p Sifthereexistsapath
i
∈ ∈{ } ∈ ∈
such that s = p and s =q. A DTMC is irreducible if any two states are reachable
1 n
fromeachother.Theperiodofastateisdefinedasthegreatestcommondivisorofthe
lengthsofallpathsthatstartandendinthestate.ADTMCisaperiodicifallstatesin
4 Mizera,PangandYuan
S are of period 1. A finite state DTMC is called ergodic if it is irreducible and aperi-
odic. By the famous ergodic theorem for DTMCs [12] an ergodic chain has a unique
stationarydistributionbeingitslimitingdistribution(alsoreferredtoasthesteady-state
distribution given by lim π Pn, where π is any initial probability distribution on
n ∞ 0 0
S).Inconsequence,thelim→itingdistributionforanergodicchainisindependentofthe
choiceofπ .Thesteady-statedistributioncanbeestimatedfromanyinitialdistribution
0
byiterativelymultiplyingitbyP.
The evolution of a first-order DTMC can be described by a stochastic recurrence
sequenceXt+1=φ(Xt,Ut+1),where Ut t N isanindependentsequenceofuniformly
distributedrealrandomvariablesover{[0,}1]∈andthetransitionfunctionφ :S [0,1] S
satisfiesthepropertythatP(φ(s,U)=s)=P foranystatess,s Sand×forany→U,
(cid:48) s,s (cid:48)
(cid:48) ∈
arealrandomvariableuniformlydistributedover[0,1].WhenSispartiallyorderedand
whenthetransitionfunctionφ(,u)ismonotonicforallu,thentheMarkovchainissaid
·
tobemonotone([13,14]).
2.2 ProbabilisticBooleannetworks(PBNs)
A PBN G(V,F) consists of a set of binary-valued nodes (also known as genes)V =
v ,v ,...,v andalistofsetsF =(F ,F ,...,F ).Foreachi 1,2,...,n theset
1 2 n 1 2 n
{ } ∈{ }
F = f(i),f(i),...,f(i) isacollectionofpredictorfunctionsfornodev,wherel(i)is
i { 1 2 l(i)} i
the number of predictor functions for v. Each f(i) F is a Boolean function defined
i j ∈ i
withrespecttoasubsetofnodesreferredtoasparentnodesofv.Thereisaprobability
i
distribution on each F F: c(i) is the probability of selecting f(i) F as the next
i ∈ j j ∈ i
predictorforvianditholdsthat∑lj(=i)1c(ji)=1.Wedenotebyvi(t)thevalueofnodeviat
timepointt N.ThestatespaceofthePBNisS= 0,1 nanditisofsize2n.Thestate
∈ { }
ofthePBNattimetisgivenbysss(t)=(v (t),v (t),...,v (t)).ThedynamicsofthePBN
1 2 n
isgivenbythesequence(sss(t))∞ .WeconsiderhereindependentPBNswherepredictor
t=0
functions for different nodes are selected independently of each other. The transition
fromsss(t)tosss(t+1)isconductedbyrandomlyselectingapredictorfunctionforeach
nodev fromF andbysynchronouslyupdatingthenodevaluesinaccordancewiththe
i i
selectedfunctions.ThereareN =∏n l(i)differentwaysinwhichthepredictorscan
i=1
be selected for all n nodes. These combinations are referred to as realisations of the
PBNandarerepresentedasn-dimensionalfunctionvectors fff =(f(1),f(2),...,f(n))
k k1 k2 kn ∈
F F ... F ,wherek 1,2,...,N andk 1,2,...,l(i) .Arealizationselected
1 2 n i
× × × ∈{ } ∈{ }
attimetisreferredtoas fff(t).Duetoindependence,P(fff )=P(fff(t)= fff )=∏n c(i).
k k i=1 ki
In PBNs with perturbations, a perturbation parameter p (0,1) is introduced to
∈
sampletheperturbationvectorγγγ(t)=(γ (t),γ (t),...,γ (t)),whereγ(t) 0,1 and
1 2 n i
P(γ(t)=1)=pforallt andi 1,2,...,n .Perturbationsprovideanalte∈rn{ative}way
i
∈{ }
toregulatethedynamicsofaPBN:thenextstateisdeterminedassss(t+1)= fff(t)(sss(t))if
γγγ(t)=000andassss(t+1)=sss(t) γγγ(t)otherwise,where istheexclusiveoroperatorfor
⊕ ⊕
vectors.Theperturbations,bythelatterupdateformula,allowthesystemtomovefrom
anystatetoanyotherstateinonesingletransition,hencerendertheunderlyingMarkov
chainirreducibleandaperiodic.Therefore,thedynamicsofaPBNwithperturbations
canbeviewedasanergodicDTMC[15].ThetransitionmatrixisgivenbyP =(1
s,s
(cid:48) −
RevivingtheTwo-stateMarkovChainApproach 5
p)n∑Nk=1111[fffk(s)=s(cid:48)]P(fffk)+(1−(1−p)n)pη(s,s(cid:48))(1−p)n−η(s,s(cid:48)),where111istheindicator
function and η(s,s) is the hamming distance between states s,s S. According to
(cid:48) (cid:48)
∈
theergodictheory,addingperturbationstoanyPBNassuresthelong-rundynamicsof
theresultingPBNisgovernedbyauniquelimitingdistribution,convergencetowhich
is independent of the choice of the initial state. However, the perturbation probability
valueshouldbechosencarefully,nottodilutethebehaviouroftheoriginalPBN.Inthis
way the ‘mathematical trick’, although introduces some noise to the original system,
allowstosignificantlysimplifytheanalysisofthesteady-statebehaviour,whichisoften
ofinterestforbiologicalsystems.
ThedensityofaPBNismeasuredwithitsfunctionnumberandparentnodesnum-
ber.ForaPBNG,itsdensityisdefinedasD(G)= 1∑NF ω(i),wherenisthenumber
n i=1
ofnodesinG,N isthetotalnumberofpredictorfunctionsinG,andω(i)isthenumber
F
ofparentnodesfortheithpredictorfunction.
Within the framework of PBNs the concept of influences is defined; it formalizes
theimpactofparentsnodesonatargetnodeandenablesitsquantification([23]).The
concept is based on the notion of a partial derivative of a Boolean function f with
respecttovariablex (1 j n):
j
≤ ≤
∂f(x)
= f(x(j,0)) f(x(j,1)),
∂xj ⊕
where isadditionmodulo2(exclusiveOR)andforl 0,1
⊕ ∈{ }
x(j,l)=(x ,x ,...,x ,l,x ,...,x ).
1 2 j 1 j+1 n
−
The influence of node x on function f is the expected value of the partial derivative
j
withrespecttotheprobabilitydistributionD(x):
(cid:104)∂f(x)(cid:105) (cid:110)∂f(x) (cid:111)
I (f)=E =P =1 =P f(x(j,0))= f(x(j,1)) .
j D
∂xj ∂xj { (cid:54) }
Letnow F bethe setof predictorsfor x withcorresponding probabilities c(i) for j=
i i j
1,...,l(i)andletI (f(i))betheinfluenceofnodex onthepredictorfunction f(i).Then,
k j k j
theinfluenceofnodex onnodex isdefinedas:
k i
l(i)
I (x)= ∑I (f(i)) c(i).
k i k j · j
j=1
Thelong-terminfluencesaretheinfluencescomputedwhenthedistributionD(x)isthe
stead-statedistributionofthePBN.
Wedefineandconsiderinthisstudytwotypesoflong-runsensitivities.
Definition1. Thelong-runsensitivitywithrespecttoselectionprobabilityperturbation
isdefinedas
s [c(i)=p]= π˜[c(i)=p] π ,
c j (cid:107) j − (cid:107)l
6 Mizera,PangandYuan
where denotesthel-norm,π isthesteady-statedistributionoftheoriginalPBN,
l
(cid:107)·(cid:107)
p [0,1] is the new value for c(i), and π˜[c(i) = p] is the steady-state probability dis-
∈ j j
tributionofthePBNperturbedasfollows.The jthselectionprobabilityfornodex is
i
replaced with c˜(i) = p and all c(i) selection probabilities for k I = 1,2,...,j
j k ∈ −j { −
1,j+1,...,l(i) arereplacedwith
}
c(i)
c˜(i)=c(i)+(c(i) p) k ,
k k j − ·∑ c(i)
l I j l
∈−
TheremainingselectionprobabilitiesoftheoriginalPBNareunchanged.
Definition2. The long-run sensitivity with respect to permanent on/off perturbations
ofanodex as
i
s [x]=max π˜[x 0] π , π˜[x 1] π ,
g i i l i l
{(cid:107) ≡ − (cid:107) (cid:107) ≡ − (cid:107)}
where π, π˜[x 0], and π˜[x 1] are the steady-state probability distributions of the
i i
originalPBN,o≡ftheoriginalP≡BNwithall f(i) F replacedby f˜(i) 0,andall f(i) F
i i
replacedby f˜(i) 1,respectively. ∈ ≡ ∈
≡
Notice that the definition of long-run sensitivity with respect to permanent on/off
perturbationsissimilarbutnotequivalenttothedefinitionoflong-runsensitivitywith
respectto1-genefunctionperturbationof[23].
3 TheTwo-stateMarkovChainApproach
3.1 Description
Werecallthetwo-stateMarkovchainapproachoriginallyintroducedin[10].Thetwo-
stateMarkovchainapproachisamethodforestimatingthesteady-stateprobabilityof
asubsetofstatesofaDTMC.InthisapproachthestatespaceofanarbitraryDTMCis
splitintotwodisjointsets,referredtoasmetastates.Oneofthemetastates,numbered
1, is the subset of interest and the other, numbered 0, is its complement. The steady-
stateprobabilityofmetastate1,denotedq,canbeestimatedbyperformingsimulations
oftheoriginalMarkovchain.Forthispurposeatwo-stateMarkovchainabstractionof
theoriginalDTMCisconsidered.Let Zt t(cid:62)0 beafamilyofbinaryrandomvariables,
{ }
where Z is the number of the meta state the original Markov chain is in at time t.
t
Zt t(cid:62)0 is a binary (0-1) stochastic process, but in general it is not a Markov chain.
{ }
However,asarguedin[10],areasonableassumptionisthatthedependencyin Zt t(cid:62)0
{ }
fallsoffrapidlywithlag.Therefore,anewprocess Zt(k) t(cid:62)0,whereZt(k)=Z1+(t 1)k,
will be approximately a first-order Markov chain fo{r k la}rge enough. A procedure−for
determiningappropriatekisgivenin[10].Thefirst-orderMarkovchainconsistsofthe
two meta states with transition probabilities α and β between them. See Figure 1 for
aconceptualillustrationoftheconstructionofthisabstraction.
Thesteady-stateprobabilityestimateqˆiscomputedfromasimulatedtrajectoryof
theoriginalDTMC.Thekeypointistodeterminetheoptimallengthofthetrajectory.
RevivingtheTwo-stateMarkovChainApproach 7
0 1
A D
�
E
1 ↵ 0 1 1 �
� �
B C ↵
(a)OriginalDTMC (b)Two-stateDTMC
Fig.1: Conceptual illustration of the idea of the two-state Markov chain construction.
(a) The state space of the original discrete-time Markov chain is split into two meta
states:statesAandBformmetastate0,whilestatesD,C,andE formmetastate1.The
splitofthestatespaceintometastatesismarkedwithdashedellipses.(b)Projectingthe
behaviouroftheoriginalchainonthetwometastatesresultsinabinary(0-1)stochastic
process.Afterpotentialsubsampling,itcanbeapproximatedasafirst-order,two-state
Markovchainwiththetransitionprobabilitiesα andβ setappropriately.
Two requirements are imposed. First, the abstraction of the DTMC, i.e., the two-state
Markov chain, should converge close to its steady-state distribution π =[π π ]. For-
0 1
mally,t satisfying P[Z(k)=i Z(k)= j] π <ε foragivenε >0andalli,j 0,1
| t | 0 − i| ∈{ }
needs to be determined.t is the so-called ‘burn-in’ period and determines the part of
the trajectory of the two-state Markov chain that needs to be discarded. Second, the
estimateqˆisrequiredtosatisfyP[q r(cid:54)qˆ(cid:54)q+r](cid:62)s,whereristherequiredpreci-
−
sionandsisaspecifiedconfidencelevel.Thisconditionisusedtodeterminethelength
of the second part of the trajectory used to compute qˆ, i.e., the sample size. Now, the
total required trajectory length of the original DTMC is then given by M+N, where
M=1+(t 1)kandN=1+( n(α,β) 1)k,wheret= m(α,β) .Thefunctionsm
− (cid:100) (cid:101)− (cid:100) (cid:101)
andndependonthetransitionsprobabilitiesα andβ andaregivenby
(cid:16) (cid:17)
m(α,β)= log mεa(xα(+αβ,β)) andn(α,β)= αβ(2−α−β)(cid:0)Φ−1(21(1+s))(cid:1)2,
log(1 α β ) (α+β)3 r2
| − − |
whereΦ 1 istheinverseofthestandardnormalcumulativedistributionfunction.The
−
expressionsformandnwereoriginallypresentedin[10].Derivationshoweverwerenot
providedandtheexpressionscontaintwooversights:intheformulaformtheabsolute
valueismissinginthedenominatorandintheformulaforntheinverseofΦshouldbe
usedinsteadofΦ.Weprovidedetailedderivationsoftheexpressionsformandninthe
AppendicesAandB,respectively.
Sinceα andβ areunknown,theyneedtobeestimated.Thisisachievediteratively
in the two-state Markov chain approach of [10]. It starts with sampling an arbitrary
initiallengthtrajectory,whichisthenusedforestimatingthevaluesofα andβ.Mand
N arecalculatedbasedontheseestimates.Next,thetrajectoryisextendedtoreachthe
required length, and α and β values are re-estimated. The new estimates are used to
8 Mizera,PangandYuan
re-calculate M and N. This process is iterated until M+N is smaller than the current
trajectory length. Finally, the resulting trajectory is used to estimate the steady-state
probabilityofmetastate1.Formoredetails,see[10].
Thetwo-stateMarkovchainapproachcanbeviewedasanaggregationmethodfor
reducingthestatespace.Generically,aggregationmethodsaggregategroupsofnodes
intheoriginalMarkovchaininaccordancewithagivenpartitionfunction,whichleads
toasmallertransitiongraph.Inprinciple,theaggregationcanbeanyMarkovchainover
this smaller transition system. What remains is the choice of the specific aggregation
whichisdeterminedbythechoiceofthetransitionprobabilities.Thepartitionfunction
is also used to obtain the so-called projection of the original process, i.e., the realisa-
tion of the original Markov chain is projected through the partition function. Ideally,
theaggregatedchainandtheprojectedrealisationshouldcoincide.However,sincethe
projection is in general not Markovian, the aggregation which is “closest” to the pro-
jection is considered instead and the “closeness” has to be appropriately defined. The
problemoffindingtheoptimalpartitionfunctioninthecasewherethedistancebetween
theprojectionandtheaggregationisquantifiedbytheKullback-Leiblerdivergencerate
(KLDR)hasbeenrecentlystudied,see,e.g.,[16,17,18].Thefocusintheseworksison
findingtheoptimalpartitionfunction,forwhichtheKLDRdistancebetweenthepro-
jectionandaggregationisminimised.Inourcasethepartitionfunctionwhichclassifies
theoriginalstatesintothetwometastatesisspecifiedbythebiologicalquestionunder
study. As shown in [17] and [18], for a given partition function, the aggregation clos-
esttotheprojectioncanbeanalyticallyobtainedprovidedthesteady-statedistribution
of the original chain is available. The steady-state probabilities are however our goal,
thus these techniques cannot be exploited for the determination of α and β transition
probabilities. The iterative, statistical estimation of the transition probabilities for the
aggregationremainstheonlyviablesolution.
3.2 Thechoiceoftheinitialsamplesize
Givengoodestimatesofα andβ,thetheoryofthetwo-stateMarkovchainpresented
aboveguaranteesthattheobtainedvaluesatisfiestheimposedprecisionrequirements.
However,thetwo-stateMarkovchainapproachstartswithgeneratingatrajectoryofthe
originalDTMCofanarbitrarilychoseninitiallength,i.e.,M +N =1+(m 1)k+
0 0 0
−
1+(n 1)k,wherem itthe‘burn-in’periodandn isthesamplesizeofthetwo-state
0 0 0
−
Markovchainabstraction.Anunfortunatechoicemayleadtofirstestimatesofα and
β that are biased and result in the new values of M and N such that M+N is either
smaller or not much larger than the initial M +N . In the former case the algorithm
0 0
stops immediately with the biased values for α, β and, more importantly, with an es-
timateforthesteady-stateprobabilitythatdoesnotsatisfytheprecisionrequirements.
Thesecondcasemayleadtothesameproblem.Asanexampleweconsideredatwo-
state Markov chain with α = 24 (0.0020214) and β = 24 (0.96). The steady-state
11873 25
probability distribution was [0.997899 0.002101]. With k =1, ε =10 6, r =10 3,
− −
s=0.95, m =5, and n =1,920 the first estimated values for α and β were 1
0 0 1918
(0.0005214) and 1, respectively. This subsequently led to M =2 and N =1,999, re-
sulting in a request for the extension of the trajectory by 76. After the extension, the
newestimatesforα andβ were 1 and1,respectively.TheseestimatesgaveM=2,
1997
RevivingtheTwo-stateMarkovChainApproach 9
N =1,920, and the algorithm stopped. The estimated steady-state probability distri-
bution was [0.99950 0.00050], which was outside the pre-specified precision interval
given byr. Independentlyrepeating theestimation for104 timesresulted in estimates
of the steady-state probabilities that were outside the pre-specified precision interval
90%oftimes.Giventheratherlargenumberofrepetitions,itcanbeconcludedthatthe
specified95%confidenceintervalwasnotreachedinthiscase.
Thereasonforthebiasedresultistheunfortunateinitialvalueforn andthefactthat
0
therealvalueofα issmall.Intheinitialisationphasethevalueofα isunderestimated
and n(α,β) calculatedbasedontheestimatedvaluesofα andβ isalmostthesame
(cid:100) (cid:101)
asn .Hence,subsequentextensionofthetrajectorydoesnotprovideanyimprovement
0
to the underestimated value of α since the elongation is short and the algorithm halts
afterthenextiteration.
Toidentifyandavoidsomeofsuchpitfalls,weconsideranumberofcasesandfor-
mulatesomeoftheconditionsinwhichthealgorithmmayfailtoachievethespecified
precision. Let n be the initial size of the sample used for initial estimation of α and
0
β. We assume that neither α nor β is zero. Then, the smallest possible estimates for
both α and β are greater than 1. Let us set an upper bound value for n to be 104.
n0 0
Formostcasesthisboundaryvalueisreasonableandweexpectn tobemuchsmaller.
0
Notice however that this is the case only if the real values of α and β are larger than
10 4.Wejustmentionherethatingeneraltheselectionofapropervalueforn heavily
− 0
depends on the real values of α and β, which are unknown a priori. From what was
statedabove,itfollowsthatbothfirstestimatesforα andβ aregreaterthan10 4.The
−
followingcasesarepossible.
– Ifbothα andβ aresmall,e.g.,lessthan0.1,thenwehavethat10 4<α,β <0.1
−
and n(α,β)>72,765 as can be seen by investigating the n(, ) function. In this
· ·
case the sample size is increased more than 7-fold which is reasonable since the
two-stateMarkovchainseemstobebad-mixingbythefirstestimatesofthevalues
forα andβ andthealgorithmasksforasignificantincreaseofthesamplesize.We
thereforeconcludethatthebad-mixingMarkovchaincasecanbeproperlyhandled
bythealgorithm.
– Both first estimates of α and β are close to 1. If α,β [0.7,0.98], the value of
∈
n(α,β)islargerthan19,000.Ifbothα,β >0.98thanthesizeofthesampledrops,
but in this case the Markov chain is highly well-mixing and short trajectories are
expectedtoprovidegoodestimates.
– The situation is somewhat different if one of the parameters is estimated to be
smallandtheotheriscloseto1asintheexampledescribedabove.Theextension
tothetrajectoryistoosmalltosignificantlychangetheestimatedvalueofthesmall
parameterandthealgorithmhalts.
Considering the above cases leads us to the observation that the following situation
needstobetreatedwithcare:Theestimatedvalueforoneoftheparametersiscloseto
1,thevalueofthesecondparameteriscloseto1,andn(α,β)iseithersmallerornot
n0
significantlylargerthann .
0
First approach: pitfall avoidance. In order to avoid this situation, we determine n
0
whichinprinciplecouldleadtoinitialinaccurateestimatesofα orβ andsuchthatthe
10 Mizera,PangandYuan
r 0.01 0.001 0.0001
s 0.9 0.95 0.975 0.9 0.95 0.975 0.9 0.95 0.975
n0 0/ [2,136] 0/ [2,1161] [2,1383] [2,1582] [2,11628] [2,13857] [2,15847]
∈
Table 1: Ranges of integer values for n that do not satisfy the ‘critical’ condition
0
n(α,β)<2n forthegivenvaluesofrands.
0
nextsamplesizegivenby n(α,β) wouldpracticallynotallowforanimprovementof
(cid:100) (cid:101)
the estimates. We reason as follows. As stated above, the ‘critical’ situation may take
place when one of the parameters is estimated to be very small, i.e., close to 1, and
n0
theincreaseinthesamplesizeisnotsignificantenoughtoimprovetheestimate.Ifthe
initialestimateisverysmall,therealvalueismostprobablyalsosmall,buttheestimate
isnotaccurate.Ifthevalueisunderestimatedtothelowestpossiblevalue,i.e., 1,on
n0
averagetheimprovementcantakeplaceonlyifthesamplesizeisincreasedatleastby
n .Therefore,withthetrade-offbetweentheaccuracyandefficiencyofthemethodin
0
mind,weproposethesamplesizetobeincreasedatleastbyn .Therefore,the‘critical’
0
situation condition is n(α,β)<2n . By analysing the function n(, ) as described in
0
· ·
detailsinAppendixD,wecandeterminethevaluesofn thatare‘safe’,i.e.,whichdo
0
notsatisfythe‘critical’condition.WepresenttheminTable1foranumberofvalues
forrands.
Second approach: controlled initial estimation of α and β. The formula for n is
asymptoticallyvalidforatwo-stateMarkovchainprovidedthatthevaluesforα andβ
areknown.However,thesevaluesarenotknownaprioriandtheyneedtobeestimated.
Unfortunately,theoriginalapproachdoesnotprovideanycontroloverthequalityofthe
initialestimateofthevaluesoftheseparameters.Incertainsituation,e.g.,asinthecase
discussedabove,thelackofsuchacontrolmechanismmayleadtoresultswithworse
statisticalconfidencelevelthanthespecifiedonegivenbys.Inthediscussedexample
s=95%,butthisvaluewasnotreachedintheperformedexperiment.Inordertoaddress
thisproblem,weproposetoextendtheinitialphaseofthetwo-stateapproachalgorithm
inthefollowingway.ThealgorithmsamplesatrajectoryofthegivenMarkovchainand
estimates the values of α and β. It might be the case that an arbitrarily chosen initial
sample size is not big enough to provide non-zero estimates for the two parameters.
If this is the case, the initial sample size is doubled and the trajectory is elongated to
collect a sample of required size. This is repeated iteratively until non-zero estimates
for both α and β are obtained. We introduce the following notation: αˆ and βˆ are the
non-zero estimates of the values of α and β, respectively. Furthermore, let n be the
0
samplesizeusedtoobtainfirstnon-zeroestimatesαˆ andβˆ,i.e.,n iseithertheinitial
0
samplesizeoritistwotopowersomemultipleoftheinitialsamplesize.
Oncethenon-zeroestimatesareavailable,thealgorithmcomputesthesamplesize
requiredtoreachthesconfidencelevelthatthetruevalueofmin(α,β)iswithinacer-
taininterval.Fordefiniteness,letusassumefromnowonthatαˆ <βˆ,whichsuggests
that min(α,β)=α. During the execution of the procedure outlined in the following
the inequality may be inverted. If this is the case, the algorithm makes corresponding
changeintheconsiderationofα andβ.