Extended G_calibrate to use the Rmpfr & gmp libaries, use the full th…

…eoretical ranges for the uniroot intervals, and work when alpha=0. Major speed-up to G_expected when alpha=0. G_variance also computed more accurately and efficiently when alpha=0. Minor speed-up to G_priorDensity for non-zero discount. Minor .version_above() and ifelse fixes. Minor speed-ups to simulation of scores/loadings (via backsolve() improvement) and speed-ups to local/column/cluster shrinkage parameters (via pre-computation). Prepared CRAN release.

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: IMIFA
 Type: Package
-Date: 2020-12-29
+Date: 2021-05-24
 Title: Infinite Mixtures of Infinite Factor Analysers and Related Models
-Version: 2.1.5
+Version: 2.1.6
 Authors@R: c(person("Keefe", "Murphy", email = "keefe.murphy@mu.ie", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-7709-3159")),
              person("Cinzia", "Viroli", email = "cinzia.viroli@unibo.it", role = "ctb", comment = c(ORCID = "0000-0002-3278-5266")),
              person("Isobel Claire", "Gormley", email = "claire.gormley@ucd.ie", role = "ctb", comment = c(ORCID = "0000-0001-7713-681X")))

R/Diagnostics.R

@@ -355,8 +355,8 @@ get_IMIFA_results.IMIFA        <- function(sims = NULL, burnin = 0L, thinning =
     GQ.temp2     <- list(AICMs = aicm, BICMs = bicm, DICs = dic)
     clust.ind    <- !any(is.element(method,   c("FA", "IFA")),
                      all(is.element(method, c("MFA", "MIFA")), G == 1))
-    sw.mx        <- ifelse(clust.ind, sw["mu.sw"],  TRUE)
-    sw.px        <- ifelse(clust.ind, sw["psi.sw"], TRUE)
+    sw.mx        <- !clust.ind || sw["mu.sw"]
+    sw.px        <- !clust.ind || sw["psi.sw"]
 
     if(!inf.Q)   {
       Q          <- if(length(n.fac)   > 1)  Q                   else n.fac

R/Gibbs_IFA.R

@@ -38,6 +38,8 @@
     Q.star       <- Q
     Q.store      <- vector("integer", n.store)
     Q.large      <- Q.big  <- FALSE
+    nu1.5        <- nu1 + 0.5
+    P.5          <- P/2
 
     mu.sigma     <- 1/sigma.mu
     mu.zero      <- as.numeric(mu.zero)
@@ -89,11 +91,11 @@
     # Shrinkage
       if(Q0) {
         load.2   <- lmat * lmat
-        phi      <- .sim_phi(Q=Q, P=P, nu1=nu1, nu2=nu2, tau=tau, load.2=load.2)
+        phi      <- .sim_phi(Q=Q, P=P, nu1.5=nu1.5, nu2=nu2, tau=tau, load.2=load.2)
         sum.term <- colSums2(phi * load.2)
         for(k in seq_len(Q)) {
-          delta[k]  <- if(k > 1) .sim_deltak(alpha.d2=alpha.d2, beta.d2=beta.d2, delta.k=delta[k], Q=Q, P=P, k=k,
-                       tau.kq=tau[k:Q], sum.term.kq=sum.term[k:Q]) else .sim_delta1(Q=Q, P=P, tau=tau, sum.term=sum.term,
+          delta[k]  <- if(k > 1) .sim_deltak(alpha.d2=alpha.d2, beta.d2=beta.d2, delta.k=delta[k], Q=Q, P.5=P.5, k=k,
+                       tau.kq=tau[k:Q], sum.term.kq=sum.term[k:Q]) else .sim_delta1(Q=Q, P.5=P.5, tau=tau, sum.term=sum.term,
                        alpha.d1=ifelse(Q1, alpha.d2, alpha.d1), beta.d1=ifelse(Q1, beta.d2, beta.d1), delta.1=delta[1L])
           tau       <- cumprod(delta)
         }

R/Gibbs_IMIFA.R

@@ -48,6 +48,8 @@
     G.store          <- vector("integer", n.store)
     act.store        <- G.store
     pi.alpha         <- cluster$pi.alpha
+    nu1.5            <- nu1 + 0.5
+    P.5              <- P/2
     if(learn.alpha) {
       alpha.store    <- ll.store
       alpha.shape    <- a.hyper[1L]
@@ -82,7 +84,7 @@
     sig.mu.sqrt      <- sqrt(sigma.mu)
     z                <- cluster$z
     nn               <- tabulate(z, nbins=trunc.G)
-    nn0              <- nn > 0
+    nn0              <- nn  > 0
     nn.ind           <- which(nn > 0)
     G.non            <- length(nn.ind)
     Q.star           <- Q
@@ -202,16 +204,16 @@
     # Shrinkage
       if(any(Q0))     {
         load.2       <- lapply(lmat[Gs], .power2)
-        phi[Gs]      <- lapply(Gs, function(g) if(n0q0[g]) .sim_phi(Q=Qs[g], P=P, nu1=nu1, nu2=nu2, tau=tau[[g]],
+        phi[Gs]      <- lapply(Gs, function(g) if(n0q0[g]) .sim_phi(Q=Qs[g], P=P, nu1.5=nu1.5, nu2=nu2, tau=tau[[g]],
                         load.2=load.2[[g]], sigma=MGPsig[g]) else .sim_phi_p(Q=Qs[g], P=P, nu1=nu1, nu2=nu2))
         sum.terms    <- lapply(Gs, function(g) if(n0q0[g]) colSums2(phi[[g]] * load.2[[g]]))
         for(g in Gs)  {
           Qg         <- Qs[g]
           Q1g        <- Q1[g]
           if(n0q0[g]) {
             for(k in seq_len(Qg)) {
-              delta[[g]][k]   <- if(k > 1) .sim_deltak(alpha.d2=alpha.d2, beta.d2=beta.d2, delta.k=delta[[g]][k], tau.kq=tau[[g]][k:Qg], P=P, Q=Qg,
-                                 k=k, sum.term.kq=sum.terms[[g]][k:Qg], sigma=MGPsig[g]) else .sim_delta1(Q=Qg, P=P, tau=tau[[g]], sum.term=sum.terms[[g]],
+              delta[[g]][k]   <- if(k > 1) .sim_deltak(alpha.d2=alpha.d2, beta.d2=beta.d2, delta.k=delta[[g]][k], tau.kq=tau[[g]][k:Qg], P.5=P.5, Q=Qg,
+                                 k=k, sum.term.kq=sum.terms[[g]][k:Qg], sigma=MGPsig[g]) else .sim_delta1(Q=Qg, P.5=P.5, tau=tau[[g]], sum.term=sum.terms[[g]],
                                  alpha.d1=ifelse(Q1g, alpha.d2, alpha.d1), beta.d1=ifelse(Q1g, beta.d2, beta.d1), delta.1=delta[[g]][1L], sigma=MGPsig[g])
               tau[[g]]        <- cumprod(delta[[g]])
             }
@@ -226,7 +228,7 @@
           nnX                 <- n0q0[Gs]
           n0Gq                <- which(nnX)
           nGq0                <- length(n0Gq)
-          MGPsig[n0Gq]        <- .sim_sigma(G=nGq0, P=P, Qs=Qs[n0Gq], rho1=rho1, rho2=rho2, sum.terms=sum.terms[n0Gq], tau=tau[n0Gq])
+          MGPsig[n0Gq]        <- .sim_sigma(G=nGq0, P.5=P.5, Qs=Qs[n0Gq], rho1=rho1, rho2=rho2, sum.terms=sum.terms[n0Gq], tau=tau[n0Gq])
           MGPsig[which(!nnX)] <- .sim_sigma_p(G=G - nGq0, rho1=rho1, rho2=rho2)
         }
       }

R/Gibbs_MIFA.R

@@ -43,6 +43,8 @@
     Q.store        <- matrix(0L, nrow=G, ncol=n.store)
     Q.large        <- Q.big <- Q.bigs <- FALSE
     err.z          <- z.err <- FALSE
+    nu1.5          <- nu1 + 0.5
+    P.5            <- P/2
 
     mu.sigma       <- 1/sigma.mu
     sig.mu.sqrt    <- sqrt(sigma.mu)
@@ -162,16 +164,16 @@
     # Shrinkage
       if(any(Q0))     {
         load.2     <- lapply(lmat, .power2)
-        phi        <- lapply(Gseq, function(g) if(n0q0[g]) .sim_phi(Q=Qs[g], P=P, nu1=nu1, nu2=nu2, tau=tau[[g]],
+        phi        <- lapply(Gseq, function(g) if(n0q0[g]) .sim_phi(Q=Qs[g], P=P, nu1.5=nu1.5, nu2=nu2, tau=tau[[g]],
                       load.2=load.2[[g]], sigma=MGPsig[g]) else .sim_phi_p(Q=Qs[g], P=P, nu1=nu1, nu2=nu2))
         sum.terms  <- lapply(Gseq, function(g) if(n0q0[g]) colSums2(phi[[g]] * load.2[[g]]))
         for(g in Gseq)  {
           Qg       <- Qs[g]
           Q1g      <- Q1[g]
           if(n0q0[g])   {
             for(k in seq_len(Qg)) {
-              delta[[g]][k] <- if(k > 1) .sim_deltak(alpha.d2=alpha.d2[g], beta.d2=beta.d2, delta.k=delta[[g]][k], tau.kq=tau[[g]][k:Qg], P=P, Q=Qg,
-                               k=k, sum.term.kq=sum.terms[[g]][k:Qg], sigma=MGPsig[g]) else .sim_delta1(Q=Qg, P=P, tau=tau[[g]], sum.term=sum.terms[[g]],
+              delta[[g]][k] <- if(k > 1) .sim_deltak(alpha.d2=alpha.d2[g], beta.d2=beta.d2, delta.k=delta[[g]][k], tau.kq=tau[[g]][k:Qg], P.5=P.5, Q=Qg,
+                               k=k, sum.term.kq=sum.terms[[g]][k:Qg], sigma=MGPsig[g]) else .sim_delta1(Q=Qg, P.5=P.5, tau=tau[[g]], sum.term=sum.terms[[g]],
                                alpha.d1=ifelse(Q1g, alpha.d2[g], alpha.d1[g]), beta.d1=ifelse(Q1g, beta.d2, beta.d1), delta.1=delta[[g]][1L], sigma=MGPsig[g])
               tau[[g]]      <- cumprod(delta[[g]])
             }
@@ -184,7 +186,7 @@
         }
         if(cluster.shrink)   {
           nGq0              <- sum(n0q0)
-          MGPsig[n0q0]      <- .sim_sigma(G=nGq0, P=P, Qs=Qs[n0q0], rho1=rho1, rho2=rho2, sum.terms=sum.terms[n0q0], tau=tau[n0q0])
+          MGPsig[n0q0]      <- .sim_sigma(G=nGq0, P.5=P.5, Qs=Qs[n0q0], rho1=rho1, rho2=rho2, sum.terms=sum.terms[n0q0], tau=tau[n0q0])
           MGPsig[!n0q0]     <- .sim_sigma_p(G=G - nGq0, rho1=rho1, rho2=rho2)
         }
       }

R/Gibbs_OMIFA.R

@@ -44,12 +44,14 @@
     Q.large          <- Q.big <- Q.bigs <- FALSE
     err.z            <- z.err <- FALSE
     G.store          <- vector("integer", n.store)
+    nu1.5            <- nu1 + 0.5
+    P.5              <- P/2
 
     mu.sigma         <- 1/sigma.mu
     sig.mu.sqrt      <- sqrt(sigma.mu)
     z                <- cluster$z
     nn               <- tabulate(z, nbins=G)
-    nn0              <- nn > 0
+    nn0              <- nn  > 0
     nn.ind           <- which(nn0)
     G.non            <- length(nn.ind)
     Q.star           <- Q
@@ -171,16 +173,16 @@
     # Shrinkage
       if(any(Q0))    {
         load.2       <- lapply(lmat, .power2)
-        phi          <- lapply(Gseq, function(g) if(n0q0[g]) .sim_phi(Q=Qs[g], P=P, nu1=nu1, nu2=nu2, tau=tau[[g]],
+        phi          <- lapply(Gseq, function(g) if(n0q0[g]) .sim_phi(Q=Qs[g], P=P, nu1.5=nu1.5, nu2=nu2, tau=tau[[g]],
                         load.2=load.2[[g]], sigma=MGPsig[g]) else .sim_phi_p(Q=Qs[g], P=P, nu1=nu1, nu2=nu2))
         sum.terms    <- lapply(Gseq, function(g) if(n0q0[g]) colSums2(phi[[g]] * load.2[[g]]))
         for(g in Gseq) {
           Qg         <- Qs[g]
           Q1g        <- Q1[g]
           if(n0q0[g])  {
             for(k in seq_len(Qg)) {
-              delta[[g]][k]   <- if(k > 1) .sim_deltak(alpha.d2=alpha.d2, beta.d2=beta.d2, delta.k=delta[[g]][k], tau.kq=tau[[g]][k:Qg], P=P, Q=Qg,
-                                 k=k, sum.term.kq=sum.terms[[g]][k:Qg], sigma=MGPsig[g]) else .sim_delta1(Q=Qg, P=P, tau=tau[[g]], sum.term=sum.terms[[g]],
+              delta[[g]][k]   <- if(k > 1) .sim_deltak(alpha.d2=alpha.d2, beta.d2=beta.d2, delta.k=delta[[g]][k], tau.kq=tau[[g]][k:Qg], P.5=P.5, Q=Qg,
+                                 k=k, sum.term.kq=sum.terms[[g]][k:Qg], sigma=MGPsig[g]) else .sim_delta1(Q=Qg, P.5=P.5, tau=tau[[g]], sum.term=sum.terms[[g]],
                                  alpha.d1=ifelse(Q1g, alpha.d2, alpha.d1), beta.d1=ifelse(Q1g, beta.d2, beta.d1), delta.1=delta[[g]][1L], sigma=MGPsig[g])
               tau[[g]]        <- cumprod(delta[[g]])
             }
@@ -193,7 +195,7 @@
         }
         if(cluster.shrink)     {
           nGq0                <- sum(n0q0)
-          MGPsig[n0q0]        <- .sim_sigma(G=nGq0, P=P, Qs=Qs[n0q0], rho1=rho1, rho2=rho2, sum.terms=sum.terms[n0q0], tau=tau[n0q0])
+          MGPsig[n0q0]        <- .sim_sigma(G=nGq0, P.5=P.5, Qs=Qs[n0q0], rho1=rho1, rho2=rho2, sum.terms=sum.terms[n0q0], tau=tau[n0q0])
           MGPsig[!n0q0]       <- .sim_sigma_p(G=G - nGq0, rho1=rho1, rho2=rho2)
         }
       }

R/IMIFA.R

@@ -6,8 +6,8 @@
 #' \itemize{
 #' \item{Type: }{Package}
 #' \item{Package: }{IMIFA}
-#' \item{Version: }{2.1.5}
-#' \item{Date: }{2020-12-29 (this version), 2017-02-02 (original release)}
+#' \item{Version: }{2.1.6}
+#' \item{Date: }{2021-05-24 (this version), 2017-02-02 (original release)}
 #' \item{Licence: }{GPL (>=2)}
 #' }
 #'

R/PlottingFunctions.R

@@ -209,7 +209,7 @@ plot.Results_IMIFA  <- function(x, plot.meth = c("all", "correlation", "density"
   if(param == "uniquenesses") {
     mat   <- switch(EXPR=uni.type, constrained=, unconstrained=mat, FALSE)
   }
-  mat     <- ifelse(n.var == 1, FALSE, mat)
+  mat     <- n.var != 1 && mat
   z.miss  <- missing(zlabels)
   if(!z.miss) {
     if(all(!is.factor(zlabels), !is.logical(zlabels), !is.numeric(zlabels)) ||
@@ -257,7 +257,7 @@ plot.Results_IMIFA  <- function(x, plot.meth = c("all", "correlation", "density"
          length(intervals) != 1))     stop("'intervals' must be a single logical indicator", call.=FALSE)
   if(any(!is.logical(mat),
          length(mat)       != 1))     stop("'mat' must be a single logical indicator",       call.=FALSE)
-  common       <- ifelse(missing(common) && all(grp.ind, !all.ind, m.sw["M.sw"], param == "scores", heat.map), FALSE, ifelse(grp.ind, common, TRUE))
+  common       <- !(missing(common) && all(grp.ind, !all.ind, m.sw["M.sw"], param == "scores", heat.map)) && (!grp.ind || common)
   if(any(!is.logical(common),
          length(common)    != 1))     stop("'common' must be a single logical indicator",    call.=FALSE)
   if(any(!is.logical(partial),
@@ -476,7 +476,7 @@ plot.Results_IMIFA  <- function(x, plot.meth = c("all", "correlation", "density"
           }
         } else   {
           base::plot(x=iter, y=x.plot[ind[1L],ind[2L],], type="l", ylab="", xlab="Iteration")
-          if(titles) graphics::title(main=list(paste0("Trace", ifelse(all.ind, ":\n", paste0(":\nLoadings - ", ifelse(grp.ind, paste0("Cluster ", g, " - "), ""))), var.names[ind[1L]], " Variable, Factor ", ind[2L])))
+          if(titles) graphics::title(main=list(paste0("Trace", ifelse(all.ind, ":\n", paste0(":\nLoadings - ",   ifelse(grp.ind, paste0("Cluster ", g, " - "), ""))), var.names[ind[1L]], " Variable, Factor ", ind[2L])))
         }
       }
 
@@ -576,7 +576,7 @@ plot.Results_IMIFA  <- function(x, plot.meth = c("all", "correlation", "density"
         } else   {
           plot.d  <- tryCatch(stats::density(x.plot[ind[1L],ind[2L],], bw="SJ"), error = function(e) stats::density(x.plot[ind[1L],ind[2L],]))
           base::plot(plot.d, main="", ylab="", xlab="")
-          if(titles) graphics::title(main=list(paste0("Density", ifelse(all.ind, ":\n", ":\nScores - "), "Observation ", obs.names[ind[1L]], ", Factor ", ind[2L])))
+          if(titles) graphics::title(main=list(paste0("Density",   ifelse(all.ind, ":\n", ":\nScores - "), "Observation ", obs.names[ind[1L]], ", Factor ", ind[2L])))
           graphics::polygon(plot.d, col=grey, border=NA)
         }
       }
@@ -604,7 +604,7 @@ plot.Results_IMIFA  <- function(x, plot.meth = c("all", "correlation", "density"
         } else   {
           plot.d  <- tryCatch(stats::density(x.plot[ind[1L],ind[2L],], bw="SJ"), error = function(e) stats::density(x.plot[ind[1L],ind[2L],]))
           base::plot(plot.d, main="", ylab="", xlab="")
-          if(titles) graphics::title(main=list(paste0("Density", ifelse(all.ind, ":\n", paste0(":\nLoadings - ", ifelse(grp.ind, paste0("Cluster ", g, " - "), ""))), var.names[ind[1L]], " Variable, Factor ", ind[2L])))
+          if(titles) graphics::title(main=list(paste0("Density",   ifelse(all.ind, ":\n", paste0(":\nLoadings - ", ifelse(grp.ind, paste0("Cluster ", g, " - "), ""))), var.names[ind[1L]], " Variable, Factor ", ind[2L])))
           graphics::polygon(plot.d, col=grey, border=NA)
         }
       }
@@ -1718,7 +1718,7 @@ plot.Results_IMIFA  <- function(x, plot.meth = c("all", "correlation", "density"
 #'
 #' #' # Other special cases of the PYP are also facilitated
 #' # G_priorDensity(N=50, alpha=c(alpha, 27.1401, 0),
-#' #                discount=c(discount, -27.1401/100, 0.8054447), type="b")
+#' #                discount=c(discount, -27.1401/100, 0.8054448), type="b")
   G_priorDensity      <- function(N, alpha, discount = 0, show.plot = TRUE, type = "h") {
     igmp       <- isNamespaceLoaded("Rmpfr")
     mpfrind    <- suppressMessages(requireNamespace("Rmpfr", quietly=TRUE)) && .version_above("gmp", "0.5-4")
@@ -1762,6 +1762,8 @@ plot.Results_IMIFA  <- function(x, plot.meth = c("all", "correlation", "density"
     if(any(discount < 0  &
           (alpha   <= 0  |
        !.IntMult(alpha, discount))))  stop("'alpha' must be a positive integer multiple of 'abs(discount)' when 'discount' is negative", call.=FALSE)
+    if(any(alpha   == 0  &
+       discount    <= 0))             stop("'discount' must be strictly positive when 'alpha'=0",   call.=FALSE)
     if(length(alpha)    != max.len) {
       alpha    <- rep(alpha,    max.len)
     }
@@ -1780,16 +1782,16 @@ plot.Results_IMIFA  <- function(x, plot.meth = c("all", "correlation", "density"
         rx[,i] <- gmp::asNumeric(abs(vnk    * Rmpfr::.bigz2mpfr(gmp::Stirling1.all(N))))
       } else    {
         if(disci > 0) {
-          vnk  <- c(Rmpfr::mpfr(0,         precBits=256),  cumsum(log(alphi + Nseq[-N]   * disci)))         -
+          vnk  <- c(Rmpfr::mpfr(0,         precBits=256),  cumsum(log(alphi + Nseq[-N]  * disci)))         -
                   log(Rmpfr::pochMpfr(alphi + 1, N - 1L)) - Nsq2 * log(disci)
         } else  {
           m    <- as.integer(alphi/abs(disci))
           mn   <- min(m, N)
           seqN <- seq_len(mn - 1L)
-          vnk  <- c(c(Rmpfr::mpfr(0,       precBits=256),  cumsum(log(m - seqN))  + seqN * log(abs(disci))) -
-                  log(Rmpfr::pochMpfr(alphi + 1, N - 1L)) - c(seqN, mn) * log(abs(disci)), rep(-Inf, N - mn))
+          vnk  <- c(c(Rmpfr::mpfr(0,       precBits=256),  cumsum(log(m - seqN)) + seqN * log(abs(disci))) -
+                  log(Rmpfr::pochMpfr(alphi + 1, N - 1L)) - c(seqN, mn) * log(abs(disci)),  rep(-Inf, N - mn))
         }
-        lnkd   <- lapply(Nseq, function(g) Rmpfr::sumBinomMpfr(g, f=function(k) Rmpfr::pochMpfr(-k * disci, N)))
+        lnkd   <- lapply(Nseq, function(g) Rmpfr::sumBinomMpfr(g, f=function(k) Rmpfr::pochMpfr(-k * disci, N), n0=1))
         rx[,i] <- gmp::asNumeric(exp(vnk    - lfactorial(Nsq2))  * abs(Rmpfr::mpfr2array(unlist(lnkd), dim=N)))
       }
     }

inst/NEWS.md

@@ -1,9 +1,14 @@
 __Infinite Mixtures of Infinite Factor Analysers__
 ==================================================  
 
+## IMIFA v2.1.6 - (_13<sup>th</sup> release [patch update]: 2021-05-24_)
 ### Bug Fixes & Miscellaneous Edits
 * Fixed breaking bugs associated with IM(I)FA slice samplers introduced in previous update.
 * `G_calibrate` function exported to augment existing `G_expected` & `G_variance` functions.
+* `G_variance` now computed more accurately and efficiently for the `alpha=0` case.
+* Major speed-up to `G_expected` for the `alpha=0` case.
+* Minor speed-ups to simulation of local/column/cluster shrinkage parameters + scores & loadings.
+* Minor speed-up to `G_priorDensity` for non-zero `discount`.
 * Minor speed-up to `psi_hyper`.
 
 ## IMIFA v2.1.5 - (_12<sup>th</sup> release [patch update]: 2020-12-29_)