Help:=proc(): print(`Sim(PrS,PrO,N), LoadedDie(L)`): print(`Seq1(K,N), ProbSS(PrO,x,s) `): print(`StupidViterbi(PrS,PrO,x) ,Viterbi1(PrS,PrO,x,s)`): print(`Viterbi(PrS,PrO,x)`): end: #Sim(PrS,PrO,N): Given a transition matrix #PrS (in the form of a list of lists) where #PrS[i][j] is the probability of switching from #state i to state j (i,j are between 1 and k:=nops(PrS)) #and PrO[i][j] means the probability of getting outcome j #(j=1.. L(:=nops(PrO[1]))) if you are currently at state i #and a positive integer N #outputs a random sequence in {1, 2, ...,L} #generated acording to this Hidden Markov Model Sim:=proc(PrS,PrO,N) local s,K,Li,i,L2: K:=nops(PrS): if not {seq(nops(PrS[i]),i=1..K)}= {K} then ERROR(`Bad input`): fi: Li:=[]: L2:=[]: s:=rand(1..K)(): for i from 1 to N do L2:=[op(L2),s]: Li:=[op(Li),LoadedDie(PrO[s])]: s:= LoadedDie(PrS[s]): od: Li,L2: end: #LoadedDie(L): throws a loaded die with nops(L) faces #given by the probability distribution L: LoadedDie:=proc(L) local m,L1,S,t,i: m:=denom(convert(L, `*`)): L1:=[seq(m*L[i],i=1..nops(L))]: S:=convert(L1,`+`): t:=rand(1..S)(): for i from 1 to nops(L1) do t:=t-L1[i]: if t<=0 then RETURN(i): fi: od: if i=nops(L1)+1 then ERROR(`Something is wrong`): fi: end: #Seq1(K,N):[1,K]^N Seq1:=proc(K,N) local i,S,j: if N=0 then RETURN({[]}): fi: S:=Seq1(K,N-1): {seq( seq( [ op(S[i]), j] , j=1..K), i=1.. nops(S) )}: end: #ProbSS(PrO,x,s): Given a transition matrix #PrS describing the Markov Process for the states #and a matrix PrO as above, a sequence x of outcomes #and a sequnce s of states, computes the probability #of getting the sequence x assuming the sequence s #Pr(x|s) ProbSS:=proc(PrO,x,s) local i: mul(PrO[s[i]][x[i]],i=1..nops(x))*1/nops(PrO): :end: #StupidViterbi(PrS,PrO,x): finds the most-likely sequence #of states that generates the sequence of outcomes x StupidViterbi:=proc(PrS,PrO,x) local champ, a,rec,S,K,N,try1: K:=nops(PrS): N:=nops(x): S:=Seq1(K,N): champ:=S[1]: rec:=ProbSS(PrO,x,champ)* mul(PrS[champ[j],champ[j+1]],j=1..nops(champ)-1): for a in S do try1:=ProbSS(PrO,x,a)*mul(PrS[a[j],a[j+1]],j=1..nops(a)-1): if try1>rec then rec:=try1: champ:=a: fi: od: champ,rec: end: #Viterbi1(PrS,PrO,x,s): Given a HMM [PrS,PrO], and #a sequence of outcomes x, and a state s #finds the champion sequence of states that is most #likely to produce x, that Ends with s Viterbi1:=proc(PrS,PrO,x,s) local K,m,champ,rec,s1,tem: option remember: K:=nops(PrS): if not (s>=1 and s<=K) then ERROR(`Bad input`): fi: if nops(x)=1 then RETURN([s],(1/K)*PrO[s][x[1]]): fi: rec:=0: for s1 from 1 to K do tem:=Viterbi1(PrS,PrO,[op(1..nops(x)-1,x)], s1)[2]*PrS[s1][s]* PrO[s][x[nops(x)]]: if tem>rec then champ:=[op(Viterbi1(PrS,PrO,[op(1..nops(x)-1,x)], s1)[1]),s]: rec:=tem: fi: od: champ,rec: end: #Viterbi(PrS,PrO,x): Given a Hidden Markov Model [PrS,PrO], and #a sequence of outcomes x #finds the champion sequence of states that is most #likely to produce x Viterbi:=proc(PrS,PrO,x) local champ,rec,tem,K,s: option remember: K:=nops(PrS): tem:=Viterbi1(PrS,PrO,x,1): champ:=tem[1]: rec:=tem[2]: for s from 2 to K do tem:=Viterbi1(PrS,PrO,x,s): if tem[2]>rec then rec:=tem[2]: champ:=tem[1]: fi: od: champ, rec: end: #Forward1(PrS,PrO,x,s): Given a HMM [PrS,PrO], and #a sequence of outcomes x, and a state s #finds the probability of x finishing at s Forward1:=proc(PrS,PrO,x,s) local K,m,s1,tem,sum1: option remember: K:=nops(PrS): if not (s>=1 and s<=K) then ERROR(`Bad input`): fi: if nops(x)=1 then RETURN((1/K)*PrO[s][x[1]]): fi: add(Forward1(PrS,PrO,[op(1..nops(x)-1,x)], s1)*PrS[s1][s]* PrO[s][x[nops(x)]],s1=1..K) end: