Comparison with depmixS4

Simple data set where depmix errors for nstates=2

x <- 1:4
model <- depmixS4::depmix(x ~ 1, nstates=2, ntimes=length(x))
depmixS4::fit(model)
#> Error in fb(init = init, A = trDens, B = dens, ntimes = ntimes(object), : NA/NaN/Inf in foreign function call (arg 10)

Data set where nstates=5 errors

data("buggy.5states", package="plotHMM")
plot(buggy.5states$logratio)

model.spec <- depmixS4::depmix(
  logratio ~ 1, data=buggy.5states, nstates=5)
set.seed(1)
model.fit <- depmixS4::fit(model.spec)
#> Error in fb(init = init, A = trDens, B = dens, ntimes = ntimes(object), : NA/NaN/Inf in foreign function call (arg 10)

Fit several data sequences

data(neuroblastoma, package="neuroblastoma")
library(data.table)
nb.dt <- data.table(neuroblastoma$profiles)
one.pro <- nb.dt[profile.id=="4" & chromosome%in%1:10]
ntimes <- rle(as.integer(one.pro$chromosome))
n.states <- 4
model.spec <- depmixS4::depmix(
  logratio ~ 1, data=one.pro,
  nstates=n.states, ntimes=ntimes$lengths)
set.seed(1)
unconstrained.fit <- depmixS4::fit(model.spec)
#> converged at iteration 155 with logLik: 1332.267
param.names <- c(mean="(Intercept)", sd="sd")
par.vec <- depmixS4::getpars(unconstrained.fit)
matrix(
  par.vec[names(par.vec) %in% param.names],
  ncol=length(param.names),
  byrow=TRUE,
  dimnames=list(state=1:n.states, parameter=names(param.names)))
#>      parameter
#> state         mean        sd
#>     1  0.022591507 0.0656462
#>     2  0.323438442 0.1177992
#>     3 -0.400759728 0.1855010
#>     4 -0.008333332 0.1170413
one.pro[, viterbi := factor(unconstrained.fit@posterior[,1]) ]
library(ggplot2)
ggplot()+
  geom_point(aes(
    position/1e6, logratio, color=viterbi),
    data=one.pro)+
  facet_grid(. ~ chromosome, scales="free", space="free")

How to constrain common variance parameter?

par.vec <- depmixS4::getpars(model.spec)
equal.groups <- rep(1, length(par.vec))
equal.groups[names(par.vec)=="sd"] <- 2
if(FALSE){
  constrained.fit <- depmixS4::fit(model.spec, equal=equal.groups)
}

Forward-backward algorithm speed comparison

one.chrom <- nb.dt[profile.id=="4" & chromosome=="2"]
n.states <- 3
model.spec <- depmixS4::depmix(
  logratio ~ 1,
  data=one.chrom,
  nstates=n.states)
log.emission.mat <- log(model.spec@dens[,1,])
log.transition.mat <- log(model.spec@trDens[1,,])
log.init.vec <- log(model.spec@init[1,])
microbenchmark::microbenchmark(depmixS4={
  result <- depmixS4::forwardbackward(model.spec)
}, plotHMM={
  fwd.list <- plotHMM::forward_interface(
    log.emission.mat, log.transition.mat, log.init.vec)
  back.mat <- plotHMM::backward_interface(
    log.emission.mat, log.transition.mat)
  mult.mat <- plotHMM::multiply_interface(
    fwd.list$log_alpha, back.mat)
  pairwise.array <- plotHMM::pairwise_interface(
    log.emission.mat, log.transition.mat,
    fwd.list$log_alpha, back.mat)
}, times=5)  
#> Unit: microseconds
#>      expr   min    lq    mean median    uq     max neval
#>  depmixS4  86.8  94.0  118.76  102.2 130.8   180.0     5
#>   plotHMM 286.8 288.2 8335.18  291.0 295.5 40514.4     5

plotHMM is 2-3x slower than depmixS4. Possibly due to (1) overhead of several function calls rather than just one, and (2) log space computations are slower than scaling.

Viterbi algorithm speed comparison

microbenchmark::microbenchmark(depmixS4={
  depmixS4::viterbi(model.spec)
}, plotHMM={
  plotHMM::viterbi_interface(
    log.emission.mat, log.transition.mat, log.init.vec)
}, times=5)
#> Unit: microseconds
#>      expr    min     lq    mean median     uq    max neval
#>  depmixS4 2732.9 4145.1 4145.80 4218.1 4758.0 4874.9     5
#>   plotHMM   11.6   13.2   31.88   13.8   29.8   91.0     5

plotHMM is about 100x faster than depmixS4, because of the overhead of loops in R (memory allocation in each iteration).