Running ‘Hasselblad1969.R’ [1s/2s]
Running the tests in ‘tests/Hasselblad1969.R’ failed.
Complete output:
  > require("SQUAREM")  # master funcion and the 4 different SQUAREM algorithms
  Loading required package: SQUAREM
  > require("setRNG")   # needed to reproduce the same results
  Loading required package: setRNG
  > 
  > # data for Poisson mixture estimation  from Hasselblad (JASA 1969)
  > poissmix.dat <- data.frame(death=0:9, freq=c(162,267,271,185,111,61,27,8,3,1))
  > 
  > y <- poissmix.dat$freq
  > tol <- 1.e-08
  > 
  > # generate a random initial guess for 3 parameters
  > setRNG(list(kind="Wichmann-Hill", normal.kind="Box-Muller", seed=123))
  > p0 <- c(runif(1),runif(2,0,4))    
  > 
  > # fixptfn = function that computes a single update of fixed-point iteration
  > # objfn = underlying objective function that needs to be maximized
  > #
  > 
  > # EM algorithm
  > 
  > poissmix.em <- function(p,y) {
  + # The fixed point mapping giving a single E and M step of the EM algorithm
  + # 
  + pnew <- rep(NA,3)
  + i <- 0:(length(y)-1)
  + zi <- p[1]*exp(-p[2])*p[2]^i / (p[1]*exp(-p[2])*p[2]^i + (1 - p[1])*exp(-p[3])*p[3]^i)
  + pnew[1] <- sum(y*zi)/sum(y)
  + pnew[2] <- sum(y*i*zi)/sum(y*zi)
  + pnew[3] <- sum(y*i*(1-zi))/sum(y*(1-zi))
  + p <- pnew
  + return(pnew)
  + }
  > 
  > poissmix.loglik <- function(p,y) {
  + # Objective function whose local minimum is a fixed point \
  + # negative log-likelihood of binary poisson mixture
  + i <- 0:(length(y)-1)
  + loglik <- y*log(p[1]*exp(-p[2])*p[2]^i/exp(lgamma(i+1)) + 
  + 		(1 - p[1])*exp(-p[3])*p[3]^i/exp(lgamma(i+1)))
  + return ( -sum(loglik) )
  + }
  > 
  > good.tol <- 1.e-05
  > 
  > good <-  c(0.6401146029910, 2.6634043566619, 1.2560951012662)
  > 
  >   pf1 <- fpiter(p=p0, y=y, fixptfn=poissmix.em, objfn=poissmix.loglik, control=list(tol=tol))
  >   print(pf1$par, digits=16)
  [1] 0.6401136228187346 2.6634055536337877 1.2560968051580201
  > #  good <- c(0.6401136228187346, 2.6634055536337877, 1.2560968051580201)
  >   if (good.tol < max(abs(good - pf1$par))) stop("error in EM algorithm.")
  > 
  > 
  > # First-order SQUAREM algorithm with SqS3 method
  >   pf2 <- squarem(par=p0, y=y, fixptfn=poissmix.em, objfn=poissmix.loglik,   
  +     control=list(tol=tol))
  >   print(pf2$par, digits=16)
  [1] 0.640115153281819 2.663403684662693 1.256094144647643
  > #[1] 0.640114475404004 2.663404513995809 1.256095320343123 linux Ubuntu 10.04 Intel core 2 duo/Dell
  > #[1] 0.640114475330512 2.663404514085158 1.256095320471585 linux Ubuntu 10.04 Intel Centrino/Toshiba
  > 
  > #  good <- c(0.640114475404004, 2.663404513995809, 1.256095320343123)
  >   if (good.tol < max(abs(good - pf2$par))) stop("error in First-order SQUAREM algorithm with SqS3 method.")
  > 
  > # First-order SQUAREM algorithm with SqS2 method
  >   pf3 <- squarem(par=p0, y=y, fixptfn=poissmix.em, objfn=poissmix.loglik, 
  +     control=list(method=2, tol=tol))
  >   print(pf3$par, digits=16)
  [1] 0.6401152282568314 2.6634035931052860 1.2560940143121484
  > #[1] 0.640115228256831  2.663403593105286  1.256094014312148 linux Ubuntu 10.04 Intel core 2 duo/Dell
  > #[1] 0.6401152283764715 2.6634035929591846 1.2560940141041679linux Ubuntu 10.04 Intel Centrino/Toshiba
  > 
  > #  good <- c(0.640115228256831, 2.663403593105286, 1.256094014312148)
  >   if (good.tol < max(abs(good - pf3$par))) stop("error in First-order SQUAREM algorithm with SqS2 method.")
  > 
  > # First-order SQUAREM algorithm with SqS3 method; non-monotone 
  > # Note: the objective function is not evaluated when objfn.inc = Inf 
  >   pf4 <- squarem(par=p0,y=y, fixptfn=poissmix.em, 
  +     control=list(tol=tol, objfn.inc=Inf))
  >   print(pf4$par, digits=16)
  [1] 0.640115153281819 2.663403684662693 1.256094144647643
  > #[1] 0.640114475404004 2.663404513995809 1.256095320343123 linux Ubuntu 10.04 Intel core 2 duo/Dell
  > #[1] 0.640114475330512 2.663404514085158 1.256095320471585 linux Ubuntu 10.04 Intel Centrino/Toshiba
  > 
  > #  good <- c(0.640114475404004, 2.663404513995809, 1.256095320343123)
  >   if (good.tol < max(abs(good - pf4$par))) stop("error in First-order SQUAREM algorithm with SqS3 method; non-monotone.")
  > 
  > # First-order SQUAREM algorithm with SqS3 method; objective function is not specified
  >   pf5 <- squarem(par=p0,y=y, fixptfn=poissmix.em, control=list(tol=tol, kr=0.1))
  >   print(pf5$par, digits=16)
  [1] 0.640115153281819 2.663403684662693 1.256094144647643
  > #[1] 0.640114475404004 2.663404513995809 1.256095320343123 linux Ubuntu 10.04 Intel core 2 duo/Dell
  > #[1] 0.640114475330512 2.663404514085158 1.256095320471585 linux Ubuntu 10.04 Intel Centrino/Toshiba
  > 
  > #  good <- c(0.640114475404004, 2.663404513995809, 1.256095320343123)
  >   if (good.tol < max(abs(good - pf5$par))) stop("error in First-order SQUAREM algorithm with SqS3 (and no objective function).")
  > 
  > # Second-order (K=2) SQUAREM algorithm with SqRRE 
  >   pf6 <- squarem(par=p0, y=y, fixptfn=poissmix.em, objfn=poissmix.loglik,
  +     control=list (K=2, tol=tol))
   ----------- FAILURE REPORT -------------- 
   --- failure: the condition has length > 1 ---
   --- srcref --- 
  : 
   --- package (from environment) --- 
  SQUAREM
   --- call from context --- 
  cyclem1(par, fixptfn, objfn, ..., control = ctrl)
   --- call from argument --- 
  if (class(p) == "try-error" | any(is.nan(unlist(p))) | any(!is.numeric(p))) stop("Error in function evaluation \n")
   --- R stacktrace ---
  where 1: cyclem1(par, fixptfn, objfn, ..., control = ctrl)
  where 2: squarem(par = p0, y = y, fixptfn = poissmix.em, objfn = poissmix.loglik, 
      control = list(K = 2, tol = tol))
  
   --- value of length: 2 type: logical ---
  [1] FALSE FALSE
   --- function from context --- 
  function (par, fixptfn, objfn, control = list(), ...) 
  {
      control.default <- list(K = 2, method = "rre", square = TRUE, 
          step.min0 = 1, step.max0 = 1, mstep = 4, kr = 1, objfn.inc = 1, 
          tol = 1e-07, maxiter = 1500, trace = FALSE, intermed = FALSE)
      namc <- names(control)
      if (!all(namc %in% names(control.default))) 
          stop("unknown names in control: ", namc[!(namc %in% names(control.default))])
      ctrl <- modifyList(control.default, control)
      method <- ctrl$method
      K <- ctrl$K
      square <- ctrl$square
      tol <- ctrl$tol
      maxiter <- ctrl$maxiter
      objfn.inc <- ctrl$objfn.inc
      trace <- ctrl$trace
      if (trace) 
          cat("Cyclem-1 \n")
      if (missing(objfn)) 
          stop("\n cyclem2 should be used if objective function is not available \n\n")
      if (K > length(par)) {
          cat("K is too large.  Decrease it.  \n")
          return()
      }
      iter <- 0
      lold <- objfn(par, ...)
      leval <- 1
      ptmp <- fixptfn(par, ...)
      res <- sqrt(c(crossprod(ptmp - par)))
      par <- ptmp
      feval <- 1
      conv <- TRUE
      while (feval < maxiter & res > tol) {
          iter <- iter + 1
          extrap <- TRUE
          p <- matrix(NA, nrow = length(par), ncol = K + 2)
          U <- matrix(NA, nrow = length(par), ncol = K + 1)
          p[, 1] <- par
          for (i in 2:(K + 2)) {
              p[, i] <- try(fixptfn(p[, i - 1], ...), silent = TRUE)
              if (class(p[, i]) == "try-error" | any(is.nan(unlist(p[, 
                  i]))) | any(!is.numeric(p[, i]))) 
                  stop("Error in function evaluation \n")
              U[, i - 1] <- p[, i] - p[, i - 1]
          }
          feval <- feval + K + 1
          sqKp1 <- sqrt(c(crossprod(U[, K + 1])))
          if (method == "rre") {
              coef <- try(solve(qr(crossprod(U), LAPACK = TRUE, 
                  tol = 1e-14), rep(1, K + 1)), silent = TRUE)
              if (class(coef) == "try-error" | any(is.nan(coef))) {
                  extrap <- FALSE
              }
              else {
                  if (abs(sum(coef)) < 1e-07) 
                    extrap <- FALSE
                  coef <- coef/sum(coef)
              }
          }
          if (method == "mpe") {
              coef <- try(solve(qr(U[, -(K + 1)], LAPACK = TRUE, 
                  tol = 1e-14), -U[, K + 1]), silent = TRUE)
              if (class(coef) == "try-error" | any(is.nan(coef))) {
                  extrap <- FALSE
              }
              else {
                  coef <- c(coef, 1)
                  if (abs(sum(coef)) < 1e-07) 
                    extrap <- FALSE
                  coef <- coef/sum(coef)
              }
          }
          if (!extrap) {
              par <- p[, ncol(p)]
              res <- sqKp1
              iter <- iter + 1
              if (trace) 
                  cat(" objective fn: ", lold, "extrapolation: ", 
                    extrap, "\n")
              next
          }
          if (!square) 
              pnew <- c(p[, -(K + 2)] %*% coef)
          if (square) {
              pnew <- rep(0, length(par))
              if (K > 1) {
                  for (i in 1:(K - 1)) {
                    p <- try(cbind(p, fixptfn(p[, i + K + 1], ...)), 
                      silent = TRUE)
                    if (class(p) == "try-error" | any(is.nan(unlist(p))) | 
                      any(!is.numeric(p))) 
                      stop("Error in function evaluation \n")
                  }
                  feval <- feval + K - 1
              }
              for (i in 0:K) {
                  for (j in 0:K) {
                    pnew <- pnew + coef[i + 1] * coef[j + 1] * 
                      p[, i + j + 1]
                  }
              }
              sqKp1 <- sqrt(c(crossprod(p[, 2 * K] - p[, 2 * K - 
                  1])))
          }
          if (extrap) {
              ptmp <- try(fixptfn(pnew, ...), silent = TRUE)
              res <- sqrt(c(crossprod(ptmp - pnew)))
              feval <- feval + 1
          }
          if (class(ptmp) == "try-error" | any(is.nan(unlist(ptmp)))) {
              pnew <- p[, ncol(p)]
              res <- sqKp1
              extrap <- FALSE
          }
          else pnew <- ptmp
          lnew <- try(objfn(pnew, ...), silent = TRUE)
          leval <- leval + 1
          if (class(lnew) == "try-error" | is.nan(lnew) | (lnew > 
              lold + objfn.inc)) {
              pnew <- p[, ncol(p)]
              res <- sqKp1
              extrap <- FALSE
              if (trace) 
                  cat(" objective fn: ", lold, "extrapolation: ", 
                    extrap, "\n")
          }
          else lold <- lnew
          par <- pnew
          if (trace) 
              cat(" objective fn: ", lold, "extrapolation: ", extrap, 
                  "\n")
      }
      if (feval >= maxiter) 
          conv <- FALSE
      return(list(par = par, value.objfn = lold, iter = iter, fpevals = feval, 
          objfevals = leval, convergence = conv))
  }
  <bytecode: 0x55e988ed0d48>
  <environment: namespace:SQUAREM>
   --- function search by body ---
  Function cyclem1 in namespace SQUAREM has this body.
   ----------- END OF FAILURE REPORT -------------- 
  Fatal error: the condition has length > 1