viterbi декодирование матриц

Question

Цель: глобальное декодирование скрытой последовательности с использованием матрицы. Код шаг за шагом ниже:

# equilibrium probs 
equil=function(P){ 
     e=eigen(t(P))$vectors[,1] 
     e/sum(e) 
} 

# simulate hidden sequence 
hssim=function(n,lambda) 
{ 
     r=dim(lambda)[1] 
     states=1:r 
     s=vector("numeric",n) 
     pi.lam=equil(lambda) 
     s[1]=sample(states,1,FALSE,pi.lam) 
     for (t in 2:n) { 
       s[t]=sample(states,1,FALSE,prob=lambda[s[t-1],]) 
     } 
     s 
} 

# simulate observed sequence 
obssim=function(s,P) 
{ 
n=length(s) 
r=dim(P)[2] 
q=dim(P)[1] 
states=1:q 
obs=vector("numeric",n) 
for (t in 1:n) { 
obs[t]=sample(states,1,FALSE,prob=P[,s[t]]) 
} 
obs 
} 

лямбда = rbind (с (0.999,0.0002,0.0003,0.0004,0.0001), с (0.001,0.9930,0.010,0.004,0.002), с (0,0004 , 0.0020,0.9900,0.0075,0.0001), с (0.0007,0.0030,0.0003,0.9940,0.0020), с (0.0010,0.0005,0.0040,0.0020,0.9925))

s=hssim(10000,lambda) 

Р = cbind (с (0,6 , 0.1,0.1,0.2), с (0.25,0.3,0.2,0.25), с (0.1,0.6,0.2,0.1), с (0.25,0.2,0.4,0.1), с (0.5,0.2,0.2,0.1))

obs=obssim(s,P) 

# optional - converting to/from another alphabet... 
letters=c("a","c","g","t") 
numbers=1:4 
convert=function(x,frm,to) 
{ 
    to[match(x,frm)] 
} 
obslets=convert(obs,numbers,letters) 

# estimate emmision probs from observed and hidden sequence 
Pest=function(s,obs,r,q) 
{ 
est=matrix(0,nrow=q,ncol=r) 
for (i in 1:r) { 
    est[,i]=table(obs[s==i]) 
    est[,i]=est[,i]/sum(est[,i]) 
    } 
    est 
} 

phat=Pest(s,obs,5,4) 

# estimate lambda from hidden sequence 
lambdaest=function(s,r) 
{ 
    n=length(s) 
    est=matrix(0,ncol=r,nrow=r) 
    for (t in 2:n) { 
     est[s[t-1],s[t]]=est[s[t-1],s[t]]+1 
     } 
    for (i in 1:r) { 
     est[i,]=est[i,]/sum(est[i,]) 
    } 
    est 
} 

lamhat=lambdaest(s,5) 
# global decoding algorithm 
global=function(obs,lambda,P) 
{ 
    r=dim(lambda)[1] 
    q=dim(P)[1] 
    n=length(obs) 
    s=vector("numeric",n) 
    f=matrix(0,nrow=r,ncol=n) 
    # forwards 
    f0=equil(lambda) 
    f[,1]=P[obs[1],]*(f0%*%lambda) 
    for (i in 2:n) { 
     for (k in 1:r){ 
     f[k,i]=P[obs[i],k]*max(f[,i-1]*lambda[,k]) 
     } 
     f[,i]=f[,i]/sum(f[,i]) 
    } 
    # backwards 
    s[n]=which.max(f[,n]) 
    for (i in (n-1):1) { 
     s[i]=which.max(lambda[,s[i+1]]*f[,i]) 
    } 
    s 
} 

globest=global(obs,lambda,P) 

Я создал функцию ион, чтобы решить viterbi-декодирование скрытой последовательности и показал ошибку, и поэтому я не мог построить график, идентифицирующий разницу. Может кто-нибудь исправить это для меня.

Answer 1

Проблема в том, что функция max не работает со сложными данными, такими как ваша матрица f. Предлагаю взять модуль перед применением макс. Пример:

x = c(1 + 1i, 2) 

# What you are doing 
max(x) 
# >> Error in max(x) : invalid 'type' (complex) of argument 

# Suggestion 
max(Mod(x)) 
# >> 2