Skip to content

Commit

Permalink
version 0.2-0
Browse files Browse the repository at this point in the history
  • Loading branch information
rolfTurner authored and gaborcsardi committed Sep 30, 2013
1 parent 37ff4cf commit bc0405f
Show file tree
Hide file tree
Showing 13 changed files with 363 additions and 225 deletions.
14 changes: 14 additions & 0 deletions ChangeLog
Expand Up @@ -254,3 +254,17 @@ Version 0.1-8 --> 0.1-9

Fixed up the dimensioning of Fortran arrays so that they
no longer use dummy 1's as the final dimension.

17 April 2013

Uploaded to CRAN

Version 0.1-9 --> 0.2-0
=======================

30 September 2013

Added the capability to use a different ispd for each
(independent) sequence in the data set.

Uploaded to CRAN
20 changes: 10 additions & 10 deletions DESCRIPTION
@@ -1,20 +1,20 @@
Package: hmm.discnp
Version: 0.1-9
Date: 2013-04-17
Version: 0.2-0
Date: 2013-09-30
Title: Hidden Markov models with discrete non-parametric observation
distributions.
Author: Rolf Turner <r.turner@auckland.ac.nz> and Limin Liu.
Maintainer: Rolf Turner <r.turner@auckland.ac.nz>
Depends: R (>= 0.99)
Description: Fits hidden Markov models with discrete non-parametric
observation distributions to data sets. Simulates data from
such models. Finds most probable underlying hidden states, the
most probable sequences of such states, and the log likelihood
of a collection of observations given the parameters of the
model.
Description: Fits hidden Markov models with discrete non-parametric
observation distributions to data sets. Simulates data
from such models. Finds most probable underlying hidden
states, the most probable sequences of such states, and the
log likelihood of a collection of observations given the
parameters of the model.
License: GPL (>= 2)
URL: http://www.math.unb.ca/~rolf/
Packaged: 2013-04-18 02:02:12 UTC; rolf
Packaged: 2013-09-30 21:51:42 UTC; rolf
NeedsCompilation: yes
Repository: CRAN
Date/Publication: 2013-04-18 07:29:04
Date/Publication: 2013-10-01 07:42:59
24 changes: 12 additions & 12 deletions MD5
@@ -1,5 +1,5 @@
2221e8fd8f99f8070a0e96dfd981edb2 *ChangeLog
1099d89012367c9dc0781cfe97decb3f *DESCRIPTION
118fb5b18c50040df12b8bd8439b432b *ChangeLog
4b191031971677ac3341e6ab481e4a3c *DESCRIPTION
e97c7b6e0d7f0eb3b8e1303451431001 *NAMESPACE
bb6f6344934b965a455a21347756e847 *R/First.R
184809b9344434a680bbc30dc96d3bb9 *R/RCS/First.R,v
Expand All @@ -22,35 +22,35 @@ f247069c6761dc5b4c228d1e33cafda8 *R/RCS/sim.hmm.R,v
90355144d471e120fc01fd1685cf445f *R/check.yval.R
ab1a44feb53f6bb9d961ef3b2d731638 *R/ffun.R
d9b85d50148cd101bed871410ac73a0d *R/fitted.hmm.discnp.R
8eb8cb180c8844a46a5d51e63cccd93c *R/hmm.R
a3382f2bb9ee6c593cd52b4016d2b52e *R/hmm.R
eb22d9be687321b5243b12879f272cee *R/init.all.R
2e3371667cb740307abd2c7b884e4cb4 *R/logLikHmm.R
b29e7ab24b3dafc8a5c8df9402fb15b9 *R/logLikHmm.R
de0b8ed3091fa489ada6c6e5355b1856 *R/mat2list.R
b92e3051eafb9511f4cae92645158ba8 *R/mps.R
0dfefa748fd129b0dbe03289e8c45f35 *R/pr.R
b3453ba428bdb359f61e253cc01e23f2 *R/recurse.R
1386236c29748ec48a8e4cc6a41649d8 *R/revise.ispd.R
3d4f088f3d59df2e9c71c09355d00545 *R/recurse.R
fced2f094ff5a74046a258e8463b590d *R/revise.ispd.R
453a3658f3432a9b1b1feda4400134d7 *R/revise.rho.R
f180f9d05c14a25f746e32f579875ff1 *R/revise.tpm.R
1343dfe5864d155602c1457e62891944 *R/sim.hmm.R
4c2d789f7ee7157a3f4c66fadbe74a7a *R/sp.R
8f4997df71f4c2f8222b1db1ae946086 *R/sp.R
dab163b5304bd23f35e2f7f9a3ec0f3a *R/viterbi.R
2656ac33803f83f5d131053efc7c4b6b *inst/READ_ME
dd41961396b612dbb4afc7fdcf73e55e *inst/Ratfor/RCS/afun.r,v
75869a842561e24d282d8ed2fc73eab8 *inst/Ratfor/RCS/bfun.r,v
f1d36db9f8e9fda32f12221e55c6dbdd *inst/Ratfor/RCS/gfun.r,v
7241424b21f65ca2f10e4dc859d31741 *inst/Ratfor/RCS/makefor,v
c7a6fbf28f918754d33b796e71d553fb *inst/Ratfor/RCS/recurse.r,v
3efff20f39b0eaa0a000afb7ffc3e277 *inst/Ratfor/RCS/recurse.r,v
770ae0f29b53d84c091de81fae008445 *inst/Ratfor/RCS/xfun.r,v
41f0e806077d282c98851bc544b9a5ea *inst/Ratfor/afun.r
17f9bbea01421df7e8c303a1fb5cbb95 *inst/Ratfor/bfun.r
5cbc32129adddcfc60272834a794ca4a *inst/Ratfor/gfun.r
40201e8652101f856140d17279cf3f4c *inst/Ratfor/makefor
d7c843cbee516d00dcc489c08cb9d069 *inst/Ratfor/recurse.r
41395d57b572e9fac85b6616aa3853e0 *inst/Ratfor/recurse.r
b97264d24cd7880ce7fbad579da6953b *inst/Ratfor/xfun.r
41aaa76a3bd7b5413aac29fe8cbb6bf2 *man/fitted.hmm.discnp.Rd
80e4f4b78b7a476a41a9f37db05a8b65 *man/hmm.Rd
3c143ea6632f17fe683cc2114326b16b *man/hmm.discnp-internal.Rd
3cbf04e92b8dad6921779aca39b73298 *man/hmm.Rd
684bf7787fac6112a443995901e90f27 *man/hmm.discnp-internal.Rd
8361273b18a9243c78bf1b8792166a51 *man/logLikHmm.Rd
f0b59556b361b40f21fb7373256de26d *man/mps.Rd
aaaaf6989cb369e176840ad9b1d8f3f9 *man/pr.Rd
Expand All @@ -60,5 +60,5 @@ d477c7563a021ac09adc74dfba50a4ca *man/sp.Rd
f1f57c6a20eb0a28ffd79efb0f8b1816 *src/afun.f
ee174ec13251743ea33ccbc6f9009291 *src/bfun.f
25ef80dc0168437aafd15bd0fd9974d8 *src/gfun.f
36f16d5a97248e7040e2fd9c59921377 *src/recurse.f
4e902426e7b897b9f3b64a295b2039e0 *src/recurse.f
dc1264652e18b677e558425ce8bd9aa2 *src/xfun.f
29 changes: 20 additions & 9 deletions R/hmm.R
@@ -1,7 +1,7 @@
hmm <- function(y,yval=NULL,par0=NULL,K=NULL,rand.start=NULL,stationary=TRUE,
mixture=FALSE,tolerance=1e-4,verbose=FALSE,itmax=200,
crit='PCLL',keep.y=TRUE,data.name=NULL)
{
hmm <- function(y,yval=NULL,par0=NULL,K=NULL,rand.start=NULL,stationary=cis,
mixture=FALSE,cis=TRUE,tolerance=1e-4,verbose=FALSE,itmax=200,
crit='PCLL',keep.y=TRUE,data.name=NULL) {
#
# Function hmm. To fit a Hidden Markov model to a data set where the
# observations come from one of a number of finite discrete
# distributions, depending on the (hidden) state of the Markov chain.
Expand All @@ -18,6 +18,12 @@ hmm <- function(y,yval=NULL,par0=NULL,K=NULL,rand.start=NULL,stationary=TRUE,
if(mixture & !stationary)
stop("Makes no sense for mixture to be non-stationary.\n")

# Check on consistancy of ``stationary'' and ``cis''.
if(stationary & !cis)
stop(paste("If the model is stationary the initial state\n",
"probability distribution must be constant\n",
"across data sequences.\n"))

# Put together a data name tag for the output.
if(is.null(data.name)) data.name <- deparse(substitute(y))

Expand Down Expand Up @@ -65,14 +71,19 @@ if(is.na(icrit)) stop('Stopping criterion not recognized.')

# Perform initial setting-up.
tpm <- par0$tpm
ispd <- revise.ispd(tpm)
if(cis) {
ispd <- revise.ispd(tpm)
} else { # Start all chains in state 1 with probability 1.
ispd <- matrix(0,K,length(y))
ispd[1,] <- 1
}
Rho <- par0$Rho
m <- nrow(Rho)
digits <- 2+ceiling(abs(log10(tolerance)))

old.theta <- c(c(tpm[,-K]),c(Rho[1:(m-1),]))
fy <- ffun(y,Rho)
rp <- recurse(fy,tpm,ispd,lns)
rp <- recurse(fy,tpm,ispd,lns,cis)
old.ll <- sum(log(rp$llc))

if(verbose) cat('\n Initial set-up completed ...\n\n')
Expand All @@ -89,7 +100,7 @@ repeat{
ispd <- if(stationary) {
revise.ispd(tpm)
} else {
revise.ispd(gamma=rp$gamma,lns=lns)
revise.ispd(gamma=rp$gamma,lns=lns,cis=cis)
}
Rho <- revise.rho(y,rp$gamma,yval)

Expand All @@ -99,7 +110,7 @@ repeat{
# to update the parameter estimates on the *next* EM
# step, if necessary).
fy <- ffun(y,Rho)
rp <- recurse(fy,tpm,ispd,lns)
rp <- recurse(fy,tpm,ispd,lns,cis)
ll <- sum(log(rp$llc))

# Test for convergence:
Expand Down Expand Up @@ -143,7 +154,7 @@ repeat{
if(length(y)==1) y <- y[[1]]
rslt <- list(Rho=Rho,tpm=tpm,ispd=ispd,log.like=ll,converged=converged,
nstep=nstep,y=if(keep.y) y else NULL, data.name=data.name,
stationary=stationary)
stationary=stationary, cis=cis)
class(rslt) <- "hmm.discnp"
rslt
}
9 changes: 6 additions & 3 deletions R/logLikHmm.R
Expand Up @@ -15,8 +15,11 @@ y <- mat2list(y)
Rho <- par$Rho
tpm <- par$tpm
ispd <- par$ispd
if(is.null(ispd)) ispd <- revise.ispd(tpm=tpm)
K <- length(ispd)
if(is.null(ispd)) {
ispd <- revise.ispd(tpm=tpm)
}
K <- length(ispd)
cis <- !(is.matrix(ispd) && ncol(ispd) > 1)

# Make sure that the entries of the vectors in y correspond
# to the row names of Rho.
Expand All @@ -27,6 +30,6 @@ if(K==1) return(sum(log(ffun(y,Rho))))

lns <- sapply(y,length)
fy <- ffun(y,Rho)
rp <- recurse(fy,tpm,ispd,lns)
rp <- recurse(fy,tpm,ispd,lns,cis)
sum(log(rp$llc))
}
5 changes: 4 additions & 1 deletion R/recurse.R
@@ -1,4 +1,4 @@
recurse <- function(fy,tpm,ispd,lns)
recurse <- function(fy,tpm,ispd,lns,cis)
{
#
# Function recurse to calculate the ``recursive probabilities'',
Expand All @@ -15,6 +15,7 @@ nxi <- L - nreps
#N <- K*M - K2
N <- K*K*nxi
epsilon <- sqrt(.Machine$double.eps)
nis <- if(cis) 1 else nreps

# Recursive probabilities:

Expand All @@ -27,6 +28,8 @@ epsilon <- sqrt(.Machine$double.eps)
epsilon=as.double(epsilon),
lns=as.integer(lns),
nstate=as.integer(K),
nis=as.integer(nis),
cis=as.logical(cis),
wrk=double(K2),
xlc=double(L),
ntot=as.integer(L),
Expand Down
64 changes: 40 additions & 24 deletions R/revise.ispd.R
@@ -1,31 +1,47 @@
revise.ispd <- function(tpm=NULL,gamma=NULL,lns=NULL) {
revise.ispd <- function(tpm=NULL,gamma=NULL,lns=NULL,cis=TRUE) {
# Function revise.ispd. To revise the initial state probability
# distribution

if(is.null(tpm) + (is.null(gamma) & is.null(lns)) != 1)
stop(paste("Either \"tpm\" should be given OR\n",
"\"gamma\" and \"lns\" should be given.\n"))

if(cis) {
if(is.null(tpm) + (is.null(gamma) & is.null(lns)) != 1)
stop(paste("When \"cis\" is TRUE either \"tpm\" should be given\n",
"OR \"gamma\" and \"lns\" should be given.\n"))

# If ispd is taken to be the steady state distribution:
if(!is.null(tpm)) {
eee <- eigen(t(tpm))
k <- match(1,round(eee$values,6))
if(length(k) != 1) {
cat('Problems with eigenvalues:\n')
print(eee$values)
stop()
}
v <- Re(eee$vectors[,k])
return(v/sum(v))
}

if(!is.null(tpm)) {
eee <- eigen(t(tpm))
k <- match(1,round(eee$values,6))
if(length(k) != 1) {
cat('Problems with eigenvalues:\n')
print(eee$values)
stop()
}
v <- Re(eee$vectors[,k])
ispd <- v/sum(v)
} else {
# Steady state not assumed:
v <- 0
j <- 1
nseq <- length(lns)
for(i in 1:nseq) {
v <- v + gamma[,j]
j <- j + lns[i]
v <- 0
j <- 1
nseq <- length(lns)
for(i in 1:nseq) {
v <- v + gamma[,j]
j <- j + lns[i]
}
ispd <- v/nseq
}
} else {
if(is.null(gamma) | is.null(lns))
stop(paste("When \"cis\" is FALSE \"gamma\" and \"lns\"\n",
"must be given.\n"))
jg <- 1
nseq <- length(lns)
K <- nrow(gamma)
ispd <- matrix(0,K,nseq)
for(j in 1:nseq) {
ispd[which.max(gamma[,jg]),j] <- 1
jg <- jg + lns[j]
}
}
v/nseq
return(ispd)
}
3 changes: 2 additions & 1 deletion R/sp.R
Expand Up @@ -6,6 +6,7 @@ sp <- function (y, object = NULL, tpm, Rho, ispd=NULL, means=FALSE)
ispd <- object$ispd
}
if(is.null(ispd)) ispd <- revise.ispd(tpm)
cis <- !(is.matrix(ispd) && ncol(ispd) > 1)
if(missing(y)) {
y <- if(!is.null(object)) object$y else NULL
if(is.null(y)) stop("No observation sequence supplied.\n")
Expand All @@ -15,7 +16,7 @@ sp <- function (y, object = NULL, tpm, Rho, ispd=NULL, means=FALSE)
if(is.null(row.names(Rho))) row.names(Rho) <- 1:nrow(Rho)
lns <- sapply(y,length)
fy <- ffun(y, Rho)
rp <- recurse(fy, tpm, ispd, lns)
rp <- recurse(fy, tpm, ispd, lns,cis)
prbs <- rp$gamma
if(means) {
yval <- ( if(!is.null(row.names(Rho)))
Expand Down

0 comments on commit bc0405f

Please sign in to comment.