首页 > 解决方案 > R中的嗯估计

问题描述

我试图通过从 2 个隐藏状态和每个状态的 3 个可能值提供的一系列观察来研究 Hmm 的转换和发射矩阵的似然估计

使用 SeqHmm 生成序列

#déclaration des matrices de transition et d 'emission

P.trans <- t(matrix(c(0.2,0.8,
                      0.7,0.3),nrow=2,ncol=2))

P.emiss <- t(matrix(c(0.2,0.6,0.2,
                      0.3,0.5,0.2),nrow=3,ncol=2))


# simuler un hmm avec la fonction du package seqHMM
sim <- simulate_hmm(n_sequence=1,transition_probs = P.trans,
                    emission_probs = P.emiss, initial_probs = c(1/2,1/2),
                    sequence_length = 100)

# transformer les observations comme un vecteur numeric 

sim.obs <- as.numeric(gsub('\\-', ' ', sim$observations))

使用 Hamilton's 1988 计算对数似然函数

llh <- function(p,y){

  # matrice de transition de la chaine cachée
  trans.mat <- matrix(c(1-p[1],p[1],
                          p[2],1-p[2]), nrow=2, ncol = 2,byrow = TRUE)

  # matrice d'émissions
  emiss.mat <- matrix(c(p[3],1-p[3]-p[4],p[4],
                          p[5], p[6],1-p[5]-p[6]), nrow=2,ncol=3,byrow = TRUE)

  # la longueur de la chaine

  T <- length(y)
  Pr <- numeric(T)

  #distribution initiale
  V <- c(1/2,1/2)

  # calcul de la vraisemblance
  for(c in 1:T){
    w = numeric(nrow(trans.mat))
    #i = Y[c-1]
    j = y[c] 
    for (k in 1:nrow(trans.mat)) {
      w[k] = 0
      for (l in 1:nrow(trans.mat)) {
        w[k] = w[k] + V[l]*trans.mat[l,k]*emiss.mat[l,j]
      }
    }
    Pr[c]= sum(w)
    V=w/sum(w)
  }

  # on ajout le - pour minimiser la fonction au lieu de la maximiser
  bjf = -sum(log(Pr))
  return(bjf)
}

使用 optim 最小化 - 对数似然函数


prob.init.optim <-c(0.5,0.1,0.3,0.4,0.3,0.7)

# lancer optim

optim(par=prob.init.optim, fn=llh, y=sim.obs) 

它返回负概率值!有人可以帮忙吗?

标签: estimation

解决方案


推荐阅读