廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

t分布の隠れマルコフモデル

今日の川柳

Rで隠れマルコフモデルをやるにはHMMpaというパッケージがあるけど、t分布は実装されていないようです。

t分布、裾が重いので外れ値に対して割と剛健です。

これが原系列です。

f:id:abrahamcow:20190204203804p:plain

色分けしてるのが未観測の状態です。

正規分布で推定した場合はこう。

f:id:abrahamcow:20190204203852p:plain

分散の大きいときと小さいときで分かれてる感じです。

t分布で推定した場合はこう。

f:id:abrahamcow:20190204203942p:plain

まあまあうまくいってます。

R のコードを貼ります。

フォワード・バックワードアルゴリズムいまいち理解しきれてなくてfb関数は勘で書いてます。

library(HMMpa)
dst <- function(x,mu,lambda,nu,log=FALSE){
  delta <- lambda*(x-mu)^2
  logdet <- log(lambda)/2
  out <- lgamma((nu+1)/2)+logdet-(log(nu*pi)/2+lgamma(nu/2)+((nu+1)/2)*log(1+delta/nu))
  if(log){
    return(out)
  }else{
    return(exp(out))
  }
}
rst <- function(n,mu,lambda,nu){
  eta <- rgamma(n,nu/2,nu/2)
  rnorm(n,mu,1/sqrt(lambda*eta))
}
softmax <- function(x){
  maxx <- max(x)
  exp(x-maxx)/sum(exp(x-maxx))
}
logsumexp <- function(x){
  maxx <- max(x)
  maxx + log(sum(exp(x-maxx)))
}
fb <- function(x, mu, lambda, nu, ln_pi, ln_A, K=2){
  El_A <- exp(ln_A)
  N <- length(x)
  lp <- matrix(0, nrow = N, ncol = K)
  for(i in 1:N){
    lp[i,] <- dst(x[i], mu, lambda, nu, log=TRUE)
  }
  ln_ZZ = array(0,dim=c(K, K, N))
  ln_alpha = matrix(0,N,K)
  ln_beta = matrix(0,N,K)
  ln_st = numeric(N)
  ln_alpha[1,] = ln_pi + lp[1,]
  ln_st[1] = logsumexp(ln_alpha[1,])
  ln_alpha[1,] = ln_alpha[1,] - ln_st[1]
  for (i in 2:N) {
    ln_alpha[i,] = as.vector(log(El_A%*%exp(ln_alpha[i-1,]))) + lp[i,]
    ln_st[i] = logsumexp(ln_alpha[i,])
    ln_alpha[i,] = (ln_alpha[i,] - ln_st[i])
  }
  for(i in (N-1):1){
    ln_beta[i,] = as.vector(log(t(El_A)%*%exp(ln_beta[i+1,] + lp[i+1,])))
    ln_beta[i,] = (ln_beta[i,] - ln_st[i+1] )
  }
  ln_Z = ln_alpha + ln_beta
  for (i in 1:(N-1)){
    ln_ZZ[,,i] = matrix(ln_alpha[i,],K,K,byrow = TRUE) + ln_A + matrix(ln_beta[i+1,] + lp[i+1,],K,K)
    ln_ZZ[,,i] = ln_ZZ[,,i] - ln_st[i+1]
  }
  return(list(s=exp(ln_Z),ss=exp(ln_ZZ),ll=sum(ln_st)))
}
EMHMMt <- function(x,Ahat,muhat,lambdahat,nuhat,maxit=1000){
  K <- length(muhat)
  dQ <- function(nu,Es,Eeta,Elogeta){
    1+sum(Es*(Elogeta-Eeta))/sum(Es)-digamma(nu/2)+log(nu/2)
  }
  N <- length(x)
  Elogpi <- log(rep(1/K,K))
  ElogA <- log(Ahat)
  ll <- numeric(maxit)
  for (j in 1:maxit) {
    tmps <- fb(x,muhat,lambdahat,nuhat,Elogpi,ElogA,K)
    Es <- tmps$s
    ll[j] <- tmps$ll
    delta <- sapply(muhat,function(mu)(x-mu)^2)
    den <- sweep(sweep(delta,2,lambdahat,"*"),2,nuhat,"+")
    Eeta <- sweep(1/den,2,nuhat+1,"*")
    Elogeta <- sweep(-log(den/2),2,digamma((nuhat+1)/2),"+")
    for (k in 1:K) {
      muhat[k] <- sum(Es[,k]*Eeta[,k]*x)/sum(Es[,k]*Eeta[,k])
      deltatmp <- (x-muhat[k])^2
      lambdahat[k] <- sum(Es[,k])/sum(Es[,k]*(deltatmp)*Eeta[,k])
      dentmp <- lambdahat*deltatmp+nuhat[k]
      Eetatmp <- (nuhat[k]+1)/dentmp
      Elogetatmp <- digamma((nuhat[k]+1)/2)-log(dentmp/2)
      solQ <-uniroot(dQ,Es=Es[,k],Eeta=Eetatmp,Elogeta=Elogetatmp,interval = c(1e-10,1e+5))
      nuhat[k] <- solQ$root 
    }
    betahat <-apply(tmps$ss,1:2,sum)
    ElogA <- sweep(log(betahat),1,log(rowSums(betahat)),"-")
  }
  EA <- sweep(betahat,1,rowSums(betahat),"/")
  list(mu=muhat,lambda=lambdahat,nu=nuhat,A=EA,s=tmps$s,ll=ll)
}

A=matrix(c(0.99,0.01,0.01,0.99),2,2,byrow = TRUE)
mu <- c(1,1.5)
lambda <- c(2,2)
nu <- c(5,5)
N <- 1000
s <- integer(N)
set.seed(1111)
s[1] <- sample.int(2,1,prob=A[1,])
for(i in 2:N){
  s[i] <- sample.int(2,1,prob=A[s[i-1],])
}
x <- rst(N,mu[s],lambda[s],nu[s])
plot(x,col=s)
ans <- EMHMMt(x,Ahat = gtools::rdirichlet(2,rep(1,2)),
              muhat=c(1,4),lambdahat =c(1,1),nuhat = c(10,10))
plot(ans$ll[-1],type="l")
ans$mu
ans$lambda
ans$nu
ans$A
plot(x,col=apply(ans$s,1,which.max))

ans_norm <- HMM_training(x,distribution_class = "norm",max_m = 2)
plot(x,col=apply(ans_norm$list_of_all_trained_HMMs[[2]]$zeta,1,which.max))