0

目的 : 行列を使用した隠れシーケンスのグローバルなデコード。ステップバイステップのコードは次のとおりです。 

# 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(c(0.999,0.0002,0.0003,0.0004,0.0001),c(0.001,0.9930,0.010,0.004,0.002),c(0.0004,0.0020,0.9900,0.0075,0.0001),c(0.0000307 ,0.9940,0.0020),c(0.0010,0.0005,0.0040,0.0020,0.9925))

s=hssim(10000,lambda)

P=cbind(c(0.6,0.1,0.1,0.2),c(0.25,0.3,0.2,0.25),c(0.1,0.6,0.2,0.1),c(0.25,0.2,0.4,0.1),c( 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)

隠しシーケンスのビタビ復号化を解決する関数を作成し、エラーを表示したため、グラフをプロットして違いを特定できませんでした.誰かが私のためにこれを修正できますか.

4

1 に答える 1

0

問題は、関数が行列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
于 2014-02-12T18:51:11.277 に答える