version 3.0-7

ChangeLog

@@ -1307,3 +1307,17 @@ when using the buggy code, but only slightly different.  And
 no crashes.  Go figure.  Have now fixed get.gl().
 
 Incremented the version number (to 3.0-6).  Uploaded to CRAN.
+
+09/02/2022
+
+Changed the response, to a decrease in the log likelihood when
+the EM algorithm is used, from throwing an error to returning the
+model corresponding to the last "valid" EM step, with a message
+pointing out that a decrease had occurred and suggesting the use
+of the returned value to provided starting values for the "bf" or
+"LM" algorithms.
+
+Adjusted the help file for hmm() accordingly.  Augmented the
+examples accordingly.
+
+Incremented the version number (to 3.0-7).  Uploaded to CRAN.

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: hmm.discnp
-Version: 3.0-6
-Date: 2020-06-27
+Version: 3.0-7
+Date: 2022-02-09
 Title: Hidden Markov Models with Discrete Non-Parametric Observation
         Distributions
 Author: Rolf Turner
@@ -19,6 +19,6 @@ LazyData: true
 ByteCompile: true
 License: GPL (>= 2)
 NeedsCompilation: yes
-Packaged: 2020-06-27 03:30:13 UTC; rolf
+Packaged: 2022-02-09 00:12:01 UTC; rolf
 Repository: CRAN
-Date/Publication: 2020-06-28 00:00:07 UTC
+Date/Publication: 2022-02-09 03:10:02 UTC

MD5

@@ -1,7 +1,7 @@
-f297d9a9da8c731ace8ed919a4ce6c01 *ChangeLog
-a0541616114d7f5c53a0bb3346331201 *DESCRIPTION
-14c3a38d432b0beee410c1255a86c006 *NAMESPACE
-af871c9d9eccb3851cf4191c01f3fcbb *R/First.R
+96c224fc05354bd93cb0e4e38f103117 *ChangeLog
+8bf709d2ca8b6c64d80395bbc0572db7 *DESCRIPTION
+e7412234af8c2326ba7d437a2467e8c3 *NAMESPACE
+c2c73afd466bfc25e1438da1592feea6 *R/First.R
 33b7623305e837aa24a6407e2484ab9f *R/RCS/First.R,v
 e098b69fcfd3119f2a7aa530ff33408a *R/RCS/check.yval.R,v
 2b70d61397d868088e2c17a567cc5f1a *R/RCS/ffun.R,v
@@ -32,20 +32,20 @@ d24a8f92658b600c1cf355b181814259 *R/cnvrtRho.R
 1e4f4753c74c51cc45a297fa7ecd0a96 *R/derivpi.R
 756d398935c0ff5d106797632605c19b *R/directMaxStuff.R
 75c2e0622b663e3d366b41abe380181d *R/extract.multipleHmmDataSets.R
-d5c3f2a60df9e26d593e6541fa324649 *R/ffun.R
+928215244277149dec1000323d8f695e *R/ffun.R
 a4fdb73a5db91e2bc024ec094d72ea3b *R/fitted.hmm.discnp.R
 8eeea906e4111fc50590cfe6eecc112b *R/fix.tpm.R
 d4203ea4eae83af910b63d38c2e4c63e *R/forgethgl.R
 b6c66fe281d103c7acac0376da26bfd0 *R/get.gl.R
 110bb9e6aea2433532f33dbaff792b94 *R/get.hgl.R
 1642fddcd75aadd3edf9624a12183610 *R/get.l.R
-8f2f6fd9955cda879436334af8d4ef17 *R/hmm.R
+d6faa7db88ae4fcb957183d9d8183f58 *R/hmm.R
 e60739a614ea30383dde5cf5cebe85f3 *R/hmmBD.R
 112ee96712a414ea4e1d7ec6915efff2 *R/hmmBI.R
-611b89ea8aec09d1451b55ea1b0f2bc4 *R/hmmLM.R
+159ada2bc32978b178080cbb2110eb82 *R/hmmLM.R
 0ee364089c031b8fe7e551ce7c902cc3 *R/hmmNumOpt.R
 eef3e3a0565fe762a93da23c875202dc *R/hmmSD.R
-3019353d6508d90879b0407be8e06525 *R/hmmUV.R
+5f42d687714a95fddb91c4a7a4c6ce06 *R/hmmUV.R
 56d938d317f960eb974ccad2de691802 *R/init.all.R
 8975e1c8b86945431c7b860967790ff3 *R/lmstep.R
 8bb98241df39bb8e577666b39c62e998 *R/logLikHmm.R
@@ -74,50 +74,50 @@ ed8d3339ab4d17f1b951068e485ad198 *R/squantCI.R
 059efc2cd64713e4f822898a593f8c66 *R/tidyList.R
 baf437c54d73a5f755c76eb6818d59ae *R/update.hmm.discnp.R
 998b9fa48aff086cdbb9e2c37e9630af *R/viterbi.R
-32954c149feb0b93b4ab0b4eec6bb632 *data/Bovine.rda
-ccf88a728ddf4c3ede54eade3f14b319 *data/Cryptosporidiosis.rda
-c0ea675d4094cce9f7f49c6a6be6a152 *data/Downloads.rda
-1a4981defc1a33fbc8d5c912b35281a6 *data/EricssonB_Jul2.rda
-535e81d6a87c437c378e8d3b8242fac8 *data/FattyLiver.rda
-dba4b56d5a063be04513a16da22046fd *data/FattyLiver2.rda
-53b08dcaca25b1ddd7f8f893119f85b5 *data/Hanta.rda
-d1a29d7c67af055cd8a063d869227431 *data/IPs.rda
-c2430277e312bab49a6ddf9ad335d78e *data/InfantEEGsleepstates.rda
-752fa81147b050a2995bebd0dc9337f7 *data/LegionnairesDisease.rda
-f9f0c661ea2fccdd1a3665f1ce169d1c *data/OffshoreRigcountsAlaska.rda
-1ee09beb7761e5d8bb56ff25356e411b *data/PriceStability.rda
-9dbbd795faaead2e1ac35eb1842cd5a1 *data/Strikes.rda
-038ee97f3598e5e30e9338f28a2e537e *data/SydColDisc.rda
-45a15243acc0646e5c7e260625fbf7a0 *data/WoodPeweeSong.rda
-8691f062ddb778278d0f4336952145d4 *data/ccprSim.rda
-153dc742a60b6a5415e282302a5bb1b2 *data/ftLiardFlows.rda
-5ec6b2939f9688e4253e5610de3feb34 *data/goldparticle380.rda
-0c0c9dee9745d81b6e308486128bf9d7 *data/lesionCount.rda
-f5afb35ca89526bdd28f51f61b80ba99 *data/linLandFlows.rda
-cc7c92cc72cf6ebf1b6b830c2ea800dc *data/portMannFlows.rda
-5e2c4f492b5ffaf43574e00b60c5391d *data/portMannSedCon.rda
-b407047f4bb669b55a3a2a1dcef44c15 *data/portMannSedLoads.rda
+b39a193f27733f56c4aab0afbfce5c36 *data/Bovine.rda
+f69a96f856090b2c0c79f53ac3a048a6 *data/Cryptosporidiosis.rda
+e633aea4d0c9aa146656350dc1fad794 *data/Downloads.rda
+59c4c6561407638c1a34be1aedbf7ce3 *data/EricssonB_Jul2.rda
+c738d7bff9378ffe95829fb222ca4d68 *data/FattyLiver.rda
+67ef1ffb0fdd629944214b36c360e7d1 *data/FattyLiver2.rda
+263a01fca207020a329f17a721463ae0 *data/Hanta.rda
+4ee1bf2ddedf4af48ba959e7b26e244c *data/IPs.rda
+dde879baa43fb44afc83675bec9facb0 *data/InfantEEGsleepstates.rda
+8f4920eabd0f1a1a6da558dc2547dbf3 *data/LegionnairesDisease.rda
+1221cc2b4d99649ad7e5dc2b55a3b596 *data/OffshoreRigcountsAlaska.rda
+85d19939e1c8ba13830a2a8df2746532 *data/PriceStability.rda
+25a599e10a3f3a3b5336b000fd196853 *data/Strikes.rda
+25228f7f449138e820fdc41cb0567bb2 *data/SydColDisc.rda
+cc48312e635dbfb64e36bf49231d6032 *data/WoodPeweeSong.rda
+9c97fafdd969db45a546d0e331fba5c3 *data/ccprSim.rda
+6164321af41bd8163dd76e28b5c91f4f *data/ftLiardFlows.rda
+4ee6732cbf67c92d357cdc6a01d45446 *data/goldparticle380.rda
+36022870b1fcc3d8994bdc54d4662559 *data/lesionCount.rda
+37cc3fcc49728feda2c1a679f88c484e *data/linLandFlows.rda
+a9484fa6529d8f1aae8752b3f1c87c3e *data/portMannFlows.rda
+54b3489e99f17532e396b104c8e49167 *data/portMannSedCon.rda
+d5d846b204bd4cc3655c392ef5426243 *data/portMannSedLoads.rda
 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
 3efff20f39b0eaa0a000afb7ffc3e277 *inst/Ratfor/RCS/recurse.r,v
 770ae0f29b53d84c091de81fae008445 *inst/Ratfor/RCS/xfun.r,v
-41f0e806077d282c98851bc544b9a5ea *inst/Ratfor/afun.r
+ee20f158a71eeafac0c2d7154397db49 *inst/Ratfor/afun.r
 07af2c190de1f7d838c3fc9a0d84e054 *inst/Ratfor/bfun.r
 21b74905cf48999c573fc16da65bd783 *inst/Ratfor/getgl.r
 92b1f7c5860d39ea7a56c24e3b9615e8 *inst/Ratfor/gethgl.r
 58fd89ffba92c2154880332f6e9c1eb2 *inst/Ratfor/getl.r
 5cbc32129adddcfc60272834a794ca4a *inst/Ratfor/gfun.r
 3ff2ab745591b58c0265a9be9114ae8f *inst/Ratfor/makefor
-7311ccb13023794fc5b8370c421e872c *inst/Ratfor/recurse.r
+16354605b2650e5a568b9a025de8034a *inst/Ratfor/recurse.r
 b97264d24cd7880ce7fbad579da6953b *inst/Ratfor/xfun.r
 ee255deea7c174b74e9c2071ab3ddd10 *man/SydColDisc.Rd
 3f39b95507c4191b458bc3f6bfc545f0 *man/anova.hmm.discnp.Rd
 2355aa790891c3e43f24eb7bc6638af1 *man/ccprSim.Rd
 0b97b2b22662765e47c4cb69a664ac3d *man/cnvrtRho.Rd
 bc86aeac058a7da2359a021b72ef5299 *man/fitted.hmm.discnp.Rd
-f0dc4cb08e21642fb5b26d32c19023af *man/hmm.Rd
+5d9be6d4b4a9dc28b093f14b10468a0d *man/hmm.Rd
 74a0394f256bb76d47c6a9eb69f96c16 *man/hmm.discnp-internal.Rd
 49b6e573c4c613dc2c33cecc0be74e73 *man/hydroDat.Rd
 380b1a58283d1143ada76b6e97f0af93 *man/lesionCount.Rd
@@ -135,12 +135,12 @@ d058d84f12f5c73cd5a2b59209275168 *man/squantCI.Rd
 896e49be313cbe6af1175c58341f80bb *man/update.hmm.discnp.Rd
 2dd71779a1ef41a4da3f76e97256539c *man/viterbi.Rd
 6e879b6b61c6f953f3a2eb8e958ac4c8 *man/weissData.Rd
-8fb5832d7d1bee2c0104215e131833c5 *src/afun.f
+948b1aff7e503c1a07620c00c1388599 *src/afun.f
 63c259f62e8ab26deb0aa8aed015b8dc *src/bfun.f
 84a6def4fee734b4891d2c817d4664d3 *src/getgl.f
 3be872b605c44cc58a8b26fd3a9b4c4b *src/gethgl.f
 f6c9de360be0cca3742c198d853437dc *src/getl.f
 d474b8baa6ca644cb43de57825558905 *src/gfun.f
 8cb448204c8c4f966e9407d708d31ff7 *src/init.c
-40f704ebb2afd93eced8725f449616ac *src/recurse.f
+f2ad3afed89d81822a82fdab6bb2aa8f *src/recurse.f
 6596740e33f4d310bc041be782455fcd *src/xfun.f

NAMESPACE

@@ -21,3 +21,4 @@ importFrom("stats", "update")
 importFrom("stats", "var")
 importFrom("stats", "quantile")
 importFrom("stats", "optimize")
+useDynLib(hmm.discnp)

R/First.R

@@ -1,7 +1,3 @@
-.onLoad <- function(lib, pkg) {
-	library.dynam("hmm.discnp", pkg, lib)
-}
-
 .onAttach <- function(lib, pkg) {
 	ver <- read.dcf(file.path(lib, pkg, "DESCRIPTION"), "Version")
         packageStartupMessage(paste(pkg, ver))

R/ffun.R

@@ -15,7 +15,8 @@ if(ncol(Y)==2) {
                                         v[is.na(v)] <- 1
                                         return(v)
                                     },yf=Y[,1],cvr=caviar)
-    return(do.call(cbind,rslt))
+    clyde <- do.call(cbind,rslt)
+    return(clyde)
 } else {
     y      <- Y$y
     Pred   <- Y[,-1,drop=FALSE]
@@ -95,11 +96,11 @@ function(Dat,Rho,type) {
 # type = 2 <--> bivariate, independent -- f2
 # type = 3 <--> bivariate, dependent   -- f3
 
-fy <- lapply(Dat,f[[type]],Rho=Rho)
-if(any(sapply(fy,is.null))) {
+fly <- lapply(Dat,f[[type]],Rho=Rho)
+if(any(sapply(fly,is.null))) {
    if(interactive) browser() else stop("Null fy component(s).\n")
 }
-fy <- do.call(rbind,fy)
+fy <- do.call(rbind,fly)
 t(fy)
 }
 })

R/hmm.R

@@ -164,7 +164,9 @@ if(K > 1) {
         }
     } else {
         levels(rslt$Rho$state) <- stnms
-        rslt <- append(rslt,after=1,list(Rho.matrix=cnvrtRho(rslt$Rho)))
+        if(identical(prednames,"Intercept")) {
+            rslt <- append(rslt,after=1,list(Rho.matrix=cnvrtRho(rslt$Rho)))
+        }
     }
 }
 

R/hmmLM.R

@@ -5,6 +5,7 @@ hmmLM <- function(y,par0,itmax=200,crit,lmc=10,tolerance,
 # hidden Markov model with discrete emissions, modelled
 # non-parametrically, using the Levenberg-Marquardt algorithm.
 #
+# discnp
 
 # Do some initial housekeeping.
 K      <- nrow(par0$tpm)

R/hmmUV.R

@@ -128,11 +128,17 @@ if(is.null(digits)) digits <- 2+ceiling(abs(log10(tolerance)))
     	chnge[1]  <- 100*(ll - old.ll)/(abs(old.ll + tolerance))
         if(ll < old.ll) {
             decll  <- round(abs(chnge[1]),digits)
-            whinge <- paste0("There was a decrease in the log likelihood in the\n",
-                             "  amount of ",decll,".  The EM algorithm\n",
-                             "  appears not to be working.  See the help for an\n",
-                             "  explanation.  Change to one of the other methods.\n")
-            stop(whinge)
+            whinge <- paste0("\n  At the ",em.step,"th step, there was a DECREASE\n",
+                             "  in the log likelihood in the amount of ",decll,".\n",
+                             "  The EM algorithm appears not to be working.\n",
+                             "  See the help for an explanation.  Change to one\n",
+                             "  of the other methods, perhaps using the current\n",
+                             "  results as starting values.\n")
+
+            message(whinge)
+            converged <- FALSE
+            nstep <- em.step-1
+            break
         }
 # Test for convergence:
         t1        <- as.vector(tpm[,-K])

inst/Ratfor/afun.r

@@ -1,4 +1,4 @@
-subroutine afun(fy,xispd,tpm,epsilon,n,nstate,wrk,xlc,alpha)
+subroutine afun(fy,xispd,tpm,epsilon,n,nstate,wrk,xlc,alpha,nr)
 implicit double precision(a-h,o-z)
 dimension wrk(nstate), xispd(nstate), xlc(n)
 dimension fy(nstate,n), tpm(nstate,nstate), alpha(nstate,n)
@@ -17,6 +17,10 @@ tsum = zero
 do j = 1,nstate {
 	wrk(j) =  fy(j,1)*xispd(j)
 	tsum = tsum + wrk(j)
+#if(nr == 1) {
+#    call dblepr("wrk(j):",-1,wrk(j),1)
+#    call dblepr("tsum:",-1,tsum,1)
+#}
 }
 
 if(tsum < epsilon) {
@@ -31,6 +35,9 @@ else {
 		alpha(j,1) = wrk(j)/tsum
 	}
 }
+kt = 1
+#call intpr("Old code; kt =",-1,kt,1)
+#call dblepr("Old code; xlc(kt)=",-1,xlc(kt),1)
 
 # Run through the remaining n-1 of the alphas (recursing!).
 do kt = 2,n {
@@ -56,6 +63,8 @@ do kt = 2,n {
 			alpha(j,kt) = wrk(j)/tsum
 		}
 	}
+#call intpr("Old code; kt =",-1,kt,1)
+#call dblepr("Old code; xlc(kt)=",-1,xlc(kt),1)
 }
 
 return

inst/Ratfor/recurse.r

@@ -25,8 +25,9 @@ do k = 1,nreps {
         }
 
 # Update the alpha's.
+nr = k
 	call afun(fy(1,kstart),xispd(1,kis),tpm,epsilon,n,nstate,wrk,
-                  xlc(kstart),alpha(1,kstart))
+                  xlc(kstart),alpha(1,kstart),nr)
 
 # Update the beta's.
 	call bfun(fy(1,kstart),tpm,epsilon,n,nstate,wrk,

man/hmm.Rd

@@ -578,29 +578,30 @@ hmm(y, yval=NULL, par0=NULL, K=NULL, rand.start=NULL,
    \code{ispd} is a named vector, the names being the same as the
    row (and column) names of \code{tpm}.
 
-   If \code{method} is equal to \code{"EM"} it \emph{may} get
-   switched to \code{"bf"} at some EM step if there is a decrease
-   in the log likelihood.  This is \dQuote{theoretically impossible}
-   but can occur in practice due to an intricacy in the way that
-   the EM algorithm treats \code{ispd} when \code{stationary}
-   is \code{TRUE}.   It turns out to be effectively impossible to
-   maximise the expected log likelihood unless the term in that
-   quantity corresponding to \code{ispd} is ignored (whence it
-   \emph{is} ignored).  Ignoring this term is \dQuote{asymptotically
-   negligible} but can have the unfortunate effect of occasionally
-   leading to a decrease in the log likelihood.
-
-   If such a decrease is detected, then the algorithm switches over
-   to using the \code{"bf"} method of maximisation (which is not
-   beset by this difficulty).  The current estimates of \code{ispd},
-   \code{tpm} and \code{Rho} are used as starting values
-   for the \code{"bf"} method.
-
-   It seems to me that it \emph{should} be the case that such
-   switching can occur only if \code{stationary} is \code{TRUE}.
-   However I have encountered instances in which the switch occurs
-   when \code{stationary} is \code{FALSE}.  I have yet to figure
-   out/track down what is going on here.
+   If \code{method} is equal to \code{"EM"} there may be a
+   \emph{decrease} (!!!) in the log likelihood at some EM step.
+   This is \dQuote{theoretically impossible} but can occur in
+   practice due to an intricacy in the way that the EM algorithm
+   treats \code{ispd} when \code{stationary} is \code{TRUE}.
+   It turns out to be effectively impossible to maximise the expected
+   log likelihood unless the term in that quantity corresponding
+   to \code{ispd} is ignored (whence it \emph{is} ignored).
+   Ignoring this term is \dQuote{asymptotically negligible} but
+   can have the unfortunate effect of occasionally leading to a
+   decrease in the log likelihood.
+
+   If such a decrease is detected, then the algorithm terminates
+   and issues a message to the effect that the decrease occurred
+   and suggesting that another method be used and that perhaps the
+   results from the current curtailed application of the \code{"EM"}
+   method could be used as starting values.
+
+   It seems to me that it \emph{should} be the case that such a
+   decrease in the log likelihood can occur only if \code{stationary}
+   is \code{TRUE}.  However I have encountered instances in which
+   a decrease occurred when \code{stationary} was \code{FALSE}.
+   I have yet to figure out/track down what is going on here.
+
 }
 \section{Note on \code{method}}{
    If the \code{method} is \code{"EM"} it is actually possible
@@ -656,11 +657,13 @@ hmm(y, yval=NULL, par0=NULL, K=NULL, rand.start=NULL,
    names default to \code{V1}, \code{V2} \dots.
 
    The probabilities of the emissions taking on their
-   various possible values are given by \deqn{\Pr(Y = y_i |
-   \boldsymbol{x}, \textrm{state}=S) = \ell_i/\sum_{j=1}^m
-   \ell_j}{\code{Pr(Y = y_i | x, state=S) = l[i]/sum(l)}} where
-   \eqn{\ell_j}{l_j} is the \eqn{j\textrm{th}}{\code{j}th} entry of
-   \eqn{\boldsymbol{\beta}^{\top}\boldsymbol{x}}{\code{t(beta)\%*\%x}}
+   various possible values are given by
+   \deqn{\Pr(Y = y_i | \boldsymbol{x}, \textrm{state}=S) =
+   \ell_i/\sum_{j=1}^m \ell_j}{%
+   \code{Pr(Y = y_i | x, state=S) = l[i]/sum(l)}}
+   where \eqn{\ell_j}{l_j} is the \eqn{j\textrm{th}}{\code{j}th}
+   entry of \eqn{\boldsymbol{\beta}^{\top}\boldsymbol{x}}{%
+   \code{t(beta)\%*\%x}}
    and  where in turn \eqn{\boldsymbol{x}}{\code{x}} is the vector
    of predictors and \eqn{\boldsymbol{\beta}}{\code{beta}} is the
    coefficient vector in the linear predicator that corresponds to
@@ -761,9 +764,21 @@ print(cnvrtRho(fit$Rho)) # A matrix of probabilities
 
 # 3.
 xxx <- with(SydColDisc,split(y,f=list(locn,depth)))
-fitSydCol <- hmm(xxx,K=2,itmax=10) # Two states: above and below the thermocline.
+Yval <- c("lo","mlo","m","mhi","hi")
+# Two states: above and below the thermocline.
+fitSydCol <- hmm(xxx,yval=Yval,K=2,verb=TRUE,itmax=10)
 
 # 4.
+\dontrun{ # Takes too long.
+fitlm <- hmm(xxx,yval=Yval,K=2,method="LM",verb=TRUE)
+fitem <- hmm(xxx,yval=Yval,K=2,verb=TRUE)
+# Algorithm terminates due to a decrease in the log likelhood
+# at EM step 64.
+newfitlm <- hmm(xxx,yval=Yval,par0=fitem,method="LM",verb=TRUE)
+# The log likelihood improves from -1900.988 to -1820.314
+}
+
+# 5.
 fitLesCount <- hmm(lesionCount,K=2,itmax=10) # Two states: relapse and remission.
 }
 \keyword{models}