version 1.1.0

cran · Feb 6, 2017 · e5bebcd · e5bebcd
1 parent ab3dbd6
commit e5bebcd
Show file tree

Hide file tree

Showing 21 changed files with 275 additions and 144 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: BAMBI
 Type: Package
 Title: Bivariate Angular Mixture Models
-Version: 1.0.1
-Date: 2017-01-03
+Version: 1.1.0
+Date: 2017-02-04
 Author: Saptarshi Chakraborty,
         Samuel W.K. Wong
 Maintainer: Saptarshi Chakraborty <c7rishi@ufl.edu>
@@ -14,6 +14,6 @@ LinkingTo: Rcpp, RcppArmadillo
 Imports: stats, graphics, Rcpp
 Depends: R (>= 3.2.0), parallel, label.switching, methods
 NeedsCompilation: yes
-Packaged: 2017-01-04 00:43:06 UTC; sap_c
+Packaged: 2017-02-06 04:55:39 UTC; sap_c
 Repository: CRAN
-Date/Publication: 2017-01-04 10:55:55
+Date/Publication: 2017-02-06 10:23:07
diff --git a/MD5 b/MD5
@@ -1,40 +1,42 @@
-9e434ce8833d8639cc5dd816b21a1a72 *DESCRIPTION
-9795ac4ccc3b9bbf9ff489246dcf536b *NAMESPACE
+89de6a25231dd3c6c350d349f3bbc9c8 *DESCRIPTION
+9446c0c22aed49892a4bc9a1adf095f1 *NAMESPACE
 b6733a9bb6815478afbd1e2baa0a01ae *R/BAMBI.R
 2b75886cb70ccc6e4c92691fbe669201 *R/RcppExports.R
 d2a4d740713be2c860b6f494a72f5231 *R/all_contour_fns.R
 b7ce7d87b46030303d30615167bb4835 *R/all_diagnostics.R
-12640d76969cba8c5c9178288fdf3f9a *R/all_model_select.R
-3ea21aa176f72da7c4a3df336978636d *R/all_postprodn_fns.R
-210481b5cf7564bf9939013473b0cbff *R/all_vm_fns.R
-f9dc74793d6cbf726198e2dbf6747c66 *R/all_vmcos_fns.R
-2ab30051cdf8a375e51c8a177bd8ec06 *R/all_vmsin_fns.R
-ff8b079f4f52808ccee59386b26a33c7 *R/all_wnorm2_fns.R
-055a90baa448a4795ec65531938f386e *R/all_wnorm_fns.R
+7f3b366cab23c4c8852bc06b3459f3cc *R/all_model_select.R
+b4830ecc7b36be612bf63d209be03608 *R/all_postprodn_fns.R
+eef49e0a3bfab8c695c3acbf376950ae *R/all_vm_fns.R
+19e2a5ca6a2d29feb8a3ae73aa9b78ea *R/all_vmcos_fns.R
+abef3e0542202ed26e23373b47cb10bf *R/all_vmsin_fns.R
+ed00ea6ee691a6359ebd13b57127c38b *R/all_wnorm2_fns.R
+67af5eabf9c8f303a0edde94d1b4b1e3 *R/all_wnorm_fns.R
 5ea045e2a3c9b94fc79723e2b4475364 *R/basic_fns.R
 4f06a73a0c2a8a9d8a48848b6ecbccd9 *R/data.R
 e974eb99f7434c59e167352b14aa380c *R/sysdata.rda
 9263837fd9658a9fb5bbdd79a64c5b3d *R/vmstart.R
 395ee0511958d872266cc44b01415b86 *R/wnormstart.R
 81772a76ae5a785d43670d78b0433f3d *data/tim8.rda
 becdb3cbeb451a1ab479993e82a305b8 *data/wind.rda
-585820303d00e7ff6c65efb5cb917d9d *man/AIC.angmcmc.Rd
+a68e2209d666680f0cb496ef364b6b31 *man/AIC.angmcmc.Rd
 d8bd3b3fba8ac9d39c84f66861db843d *man/DIC.Rd
 cb0fafc768e53ce379e353076598c4b8 *man/WAIC.Rd
+3f8c5ec02b9e20797cf7bf9cb50cc13a *man/bestmodel.Rd
 aebc9cc74da0f47d324fed648e71eb3a *man/contour.angmcmc.Rd
-e7508df5a72a7eb7a60dfe20445eb744 *man/d_fitted.Rd
+d4bb88a32bcaea258db1e6554fb72d80 *man/d_fitted.Rd
 c9a4dd46524559d762ab04800962c7a0 *man/densityplot1d.Rd
 ea4e49e1ccec872a4f45eeeae9d9c31c *man/densityplot2d.Rd
 925f57f5fac6444ed6b98004e763f193 *man/extractsamples.Rd
-3c54189f7f356d83a822a821b6361854 *man/fit_stepwise_bivariate.Rd
-f580b2b2805638d803df67616a6f2b34 *man/fit_stepwise_univariate.Rd
-cbb5475c030a57c982d9d6a3ede982b5 *man/fit_vmcosmix.Rd
-470f6bdfc7acd71192cb2ac9a852200c *man/fit_vmmix.Rd
-cb7b93bdc455e16ff318fc33c9f143ce *man/fit_vmsinmix.Rd
-5226e211945ae6fb04e45c04bd45bdf1 *man/fit_wnorm2mix.Rd
-6ed2d2bf576dc2c710baf241a4e22819 *man/fit_wnormmix.Rd
+5a3848acd6c77a5edf64bdcd3b4c5fcd *man/fit_stepwise_bivariate.Rd
+8b7b05b0d38fbec555fa045dcf6cc730 *man/fit_stepwise_univariate.Rd
+822aea5db0b7cb94a3d4eb9b1dc950a6 *man/fit_vmcosmix.Rd
+87280f0c9730de775e4b50e275b6c42d *man/fit_vmmix.Rd
+7fbbfa2dd9476eea0dfc1c2557439642 *man/fit_vmsinmix.Rd
+67216a4c685bf28e2d33cd0af88db439 *man/fit_wnorm2mix.Rd
+1b5516ea85fedae7310418d642ed7a05 *man/fit_wnormmix.Rd
 8b452486af9d6c261622807df43ab543 *man/fix_label.Rd
 0d10af677de2d5f18f8d0738ac8332fe *man/is.angmcmc.Rd
+42c5da279b5bf3a0d684490afdc43f09 *man/logLik.angmcmc.Rd
 042a2dfb786a28dc47c30dbe28a1a639 *man/lpdtrace.Rd
 766ba2fda125f633c9adc39c83aef940 *man/paramtrace.Rd
 397dc620107fc800e7a6a514cd78d44c *man/pointest.Rd

diff --git a/NAMESPACE b/NAMESPACE
@@ -3,12 +3,14 @@
 S3method(AIC,angmcmc)
 S3method(BIC,angmcmc)
 S3method(contour,angmcmc)
+S3method(logLik,angmcmc)
 S3method(print,angmcmc)
 S3method(print,stepfit)
 S3method(quantile,angmcmc)
 S3method(summary,angmcmc)
 export(DIC)
 export(WAIC)
+export(bestmodel)
 export(d_fitted)
 export(densityplot1d)
 export(densityplot2d)

diff --git a/R/all_model_select.R b/R/all_model_select.R
@@ -43,7 +43,7 @@
 #' fit.vm.step.15 <- fit_stepwise_univariate(wind, "vm", start_ncomp = 1,
 #'                                           max_ncomp = 3, n.iter = 15,
 #'                                           ncores = 1)
-#' (fit.vm.best.15 <- fit.vm.step.15$fit.best)
+#' (fit.vm.best.15 <- bestmodel(fit.vm.step.15))
 #' densityplot1d(fit.vm.best.15)
 #'
 #' @export
@@ -150,7 +150,7 @@ fit_stepwise_univariate <- function(data, model, fn = mean,  start_ncomp=1, max_
 #' fit.vmsin.step.15 <- fit_stepwise_bivariate(tim8, "vmsin", start_ncomp = 3,
 #'                                             max_ncomp = 5, n.iter = 15,
 #'                                             ncores = 1)
-#' (fit.vmsin.best.15 <- fit.vmsin.step.15$fit.best)
+#' (fit.vmsin.best.15 <- bestmodel(fit.vmsin.step.15))
 #' contour(fit.vmsin.best.15)
 #'
 #' @export
@@ -214,3 +214,26 @@ fit_stepwise_bivariate <- function(data, model, fn = mean, start_ncomp = 1, max_
 
   result
 }
+
+
+#' Extracting angmcmc object corresponding to the best fitted model in stepwise fits
+#'
+#' @param step_object stepwise fitted object (output of \code{\link{fit_stepwise_univariate}}
+#' or \code{\link{fit_stepwise_bivariate}}).
+#'
+#' @return Returns an angmcmc object corresponding to the the best fitted model in step_object.
+#'
+#' @examples
+#' # illustration only - more iterations needed for convergence
+#' fit.vmsin.step.15 <- fit_stepwise_bivariate(tim8, "vmsin", start_ncomp = 3,
+#'                                             max_ncomp = 5, n.iter = 15,
+#'                                             ncores = 1)
+#' fit.vmsin.best.15 <- bestmodel(fit.vmsin.step.15)
+#' fit.vmsin.best.15
+#'
+#' @export
+
+bestmodel <- function(step_object) {
+  if(class(step_object) != "stepfit") stop("\'step_object\' is not a stepwise fitted object")
+  step_object$fit.best
+}
diff --git a/R/all_postprodn_fns.R b/R/all_postprodn_fns.R
@@ -285,7 +285,7 @@ summary.angmcmc <- function(object, burnin = 1/3, thin = 1, ...)
 #' @details
 #' Let \eqn{\hat{L}} be the maximum value of the likelihood function for the model, \eqn{m} be the number of
 #' estimated parameters in the model and \eqn{n} be the number of data points. Then AIC and BIC are defined as
-#' \eqn{AIC = -2 \log \hat{L} + 2m} and \eqn{BIC = -2 \log \hat{L} + m \log(n)}.
+#' \eqn{AIC = -2 \log \hat{L} + mk} and \eqn{BIC = -2 \log \hat{L} + m \log(n)}.
 #'
 #' \eqn{\hat{L}} is estimated by the sample maximum obtained from the MCMC realizations.
 #'
@@ -298,7 +298,7 @@ summary.angmcmc <- function(object, burnin = 1/3, thin = 1, ...)
 #'
 #' @export
 
-AIC.angmcmc <- function(object, burnin = 1/3, thin = 1, ...)
+AIC.angmcmc <- function(object, burnin = 1/3, thin = 1, k = 2, ...)
 {
   if(!is.null(object$fixed.label)) {
     if(missing(burnin)) burnin <- object$burnin
@@ -314,7 +314,7 @@ AIC.angmcmc <- function(object, burnin = 1/3, thin = 1, ...)
 
   llik <- max(object$llik[final_iter])
   npar <- prod(dim(object$par.value)[1:2]) - 1
-  aic <- 2 * (npar - llik)
+  aic <- k * npar - 2 * llik
   aic
 }
 
@@ -516,14 +516,14 @@ WAIC <- function(object, form = 1, burnin = 1/3, thin = 1, ...)
 
     if(ncomp > 1) {
       all_den_mat <- (sapply(1:length(final_iter),
-                                             function(iter) do.call(paste0("d", object$model, "mix"),
-                                                                    addtolist(list_by_row(all_par[ , , iter], 1:ncomp),
-                                                                              x = object$data, int.displ = int.displ) )))
+                             function(iter) do.call(paste0("d", object$model, "mix"),
+                                                    addtolist(list_by_row(all_par[ , , iter], 1:ncomp),
+                                                              x = object$data, int.displ = int.displ) )))
     } else {
       all_den_mat <- (sapply(1:length(final_iter),
-                                             function(iter) do.call(paste0("d", object$model),
-                                                                    addtolist(as.list(all_par[-1 , iter]),
-                                                                              x = object$data, int.displ = int.displ) )))
+                             function(iter) do.call(paste0("d", object$model),
+                                                    addtolist(as.list(all_par[-1 , iter]),
+                                                              x = object$data, int.displ = int.displ) )))
     }
 
   }
@@ -555,7 +555,7 @@ WAIC <- function(object, form = 1, burnin = 1/3, thin = 1, ...)
 #'
 #' @details
 #' To estimate the mixture density, first the parameter vector \eqn{\eta} is estimated
-#' by applying \code{fn} on the MCMC samples, yielding the (consistent) Bayes estimate \eqn{\hat{\eta}}. Then the mixture density
+#' by applying \code{fn} on the MCMC samples (using the function \link{pointest}), yielding the (consistent) Bayes estimate \eqn{\hat{\eta}}. Then the mixture density
 #' \eqn{f(x|\eta)} at any point \eqn{x} is (consistently) estimated by \eqn{f(x|\hat{\eta})}.
 #'
 #' The random deviates are generated from the estimated mixture density \eqn{f(x|\hat{\eta})}.
@@ -574,9 +574,10 @@ WAIC <- function(object, form = 1, burnin = 1/3, thin = 1, ...)
 d_fitted <- function(x, object, fn = mean, burnin = 1/3, thin = 1)
 {
   if(class(object) != "angmcmc") stop("object must be an angmcmc object")
-  if((length(dim(x)) < 2 && length(x) != 2) || (length(dim(x)) == 2 && tail(dim(x), 1) != 2)
-     || (length(dim(x)) > 2)) stop("x must either be a bivariate vector or a two-column matrix")
-
+  if(object$type == "bi") {
+    if((length(dim(x)) < 2 && length(x) != 2) || (length(dim(x)) == 2 && tail(dim(x), 1) != 2)
+       || (length(dim(x)) > 2)) stop("x must either be a bivariate vector or a two-column matrix")
+  }
   if(missing(burnin)) {
     if(is.null(object$fixed.label)) burnin <- 1/3
     else burnin <- object$burnin
@@ -647,3 +648,51 @@ r_fitted <- function(n, object, fn = mean, burnin = 1/3, thin = 1)
 
 }
 
+
+#' Extract Log-Likelihood from angmcmc objects
+#' @inheritParams pointest
+#' @param fn function to evaluate on MCMC samples to estimate parameters.  Defaults to \code{mean}, which computes the estimated posterior mean. Used for parameter estimatation.
+#' @param burnin initial fraction of the MCMC samples to be discarded as burn-in. Must be a value in [0, 1). Used for parameter estimatation.
+#' @param thin positive integer. If \code{thin =} \eqn{n}, only every \eqn{n}-th realizations of the Markov chain is kept. Used for parameter estimatation.
+#' @inheritParams stats::logLik
+#'
+#' @details In order to estimate the log likelihood for the model, first the parameter vector is estimated using \link{pointest},
+#' and then log the likelihood is calculated on the basis of the estimated parameter.
+#'
+#' The degrees of the likelihood function is the total number of free parameters estimated in the mixture models,
+#' which is equal to \eqn{6K - 1} for bivariate models (vmsin, vmcos and wnorm2), or \eqn{3K - 1} for univariate
+#' models (vm and wnorm), where \eqn{K} denotes the number of components in the mixture model.
+#'
+#' @return Returns an object of class \link{logLik}. This is a number (the estimated log likelihood) with attributes "df"
+#' (degrees of freedom) and "nobs" (number of observations).
+#'
+#' @examples
+#' # illustration only - more iterations needed for convergence
+#' fit.vmsin.20 <- fit_vmsinmix(tim8, ncomp = 3, n.iter =  20,
+#'                              ncores = 1)
+#' logLik(fit.vmsin.20)
+#' @export
+
+
+logLik.angmcmc <- function(object, fn = mean, burnin = 1/3, thin = 1, ...)
+{
+  if(class(object) != "angmcmc") stop("\"object\" must be an angmcmc object")
+  if(!is.null(object$fixed.label)) {
+    if(missing(burnin)) burnin <- object$burnin
+    if(missing(thin)) thin <- object$thin
+  }
+  if(burnin < 0 || burnin >= 1) stop("\"burnin\" must be in [0, 1)")
+  if(thin <= 0  || thin != as.integer(thin) ) stop("\"thin\" must be a positive integer")
+
+  all_den <- d_fitted(object$data, object, fn = fn, burnin = burnin, thin = thin)
+  llik <- sum(log(all_den))
+
+  if(object$type == "uni") df <- 3*object$ncomp - 1
+  else df <- 6*object$ncomp - 1
+
+  object_nobs <- object$n.data
+  result <- llik
+  attributes(result) <- list("df" = df, "nobs" = object_nobs)
+  class(result) <- "logLik"
+  result
+}
diff --git a/R/all_vm_fns.R b/R/all_vm_fns.R
@@ -200,13 +200,6 @@ dvmmix <- function(x, kappa, mu, pmix)
 #' @param gam.loc,gam.scale location and scale (hyper-) parameters for the gamma prior for \code{kappa}. See
 #' \link{dgamma}. Defaults are \code{gam.loc = 0, gam.scale = 1000} that makes the prior non-informative.
 #'
-#' @usage
-#' fit_vmmix(data, ncomp, start_par = list(), method = "hmc",
-#'           epsilon = 0.07, L = 10, epsilon.random = TRUE,
-#'           L.random = FALSE, propscale = rep(0.01, 2),
-#'           n.iter = 500, gam.loc = 0, gam.scale = 1000,
-#'           autotune = FALSE, iter.tune=10, ncores,
-#'           show.progress = TRUE)
 #' @return returns an angular MCMC object.
 #'
 #' @details
@@ -228,7 +221,7 @@ dvmmix <- function(x, kappa, mu, pmix)
 
 
 fit_vmmix <- function(data, ncomp, start_par = list(), method = "hmc", epsilon=0.07, L=10, epsilon.random=TRUE,
-                      L.random=FALSE, propscale = rep(0.01, 2), n.iter=500, gam.loc=0, gam.scale=1000,
+                      L.random=FALSE, propscale = rep(0.01, 2), n.iter=500, gam.loc=0, gam.scale=1000, pmix.alpha = 1/2,
                       autotune = FALSE, iter.tune=10, ncores, show.progress = TRUE) {
 
   if(!(mode(data) %in% c("numeric", "list"))) stop("non-compatible data")
@@ -279,7 +272,7 @@ fit_vmmix <- function(data, ncomp, start_par = list(), method = "hmc", epsilon=0
   starting$par.mat[abs(starting$par.mat) >= kappa_upper/2] <- kappa_upper/2
   starting$l.c.univm <- as.numeric(log_const_univm_all(starting$par.mat))
   starting$llik <- llik_univm_full(data.rad, starting$par.mat, starting$pi.mix, starting$l.c.univm, ncores)
-  starting$lprior <- sum(ldgamanum(starting$par.mat[1,], 0, 1000))
+  starting$lprior <- sum((pmix.alpha-1) * log(starting$pi.mix)) + sum(ldgamanum(starting$par.mat[1,], gam.loc, gam.scale))
   starting$lpd <- starting$llik + starting$lprior
 
   par.mat.all <- array(0, dim = c(2, ncomp, n.iter+1))
@@ -465,7 +458,7 @@ fit_vmmix <- function(data, ncomp, start_par = list(), method = "hmc", epsilon=0
         par.mat.prop <- q
         l.c.univm.prop <- log_const_univm_all(par.mat.prop)
 
-        lprior.prop <- sum(ldgamanum(q[1,], gam.loc, gam.scale))
+        lprior.prop <- sum((pmix.alpha-1) * log(pi.mix.1)) + sum(ldgamanum(q[1,], gam.loc, gam.scale))
 
         llik.prop <- llik_univm_full(data.rad, q, pi.mix.1, l.c.univm.prop, ncores)
 
@@ -541,7 +534,7 @@ fit_vmmix <- function(data, ncomp, start_par = list(), method = "hmc", epsilon=0
       lprior_old <- MC$lprior
 
       llik_prop <- llik_univm_full(data.rad, prop.mat, pi.mix.1, l.c.univm.prop, ncores)
-      lprior_prop <- sum(ldgamanum(k.1.prop, 0, 1000))
+      lprior_prop <- sum((pmix.alpha-1) * log(pi.mix.1)) + sum(ldgamanum(k.1.prop, gam.loc, gam.scale))
 
       lpd_old <- llik_old + lprior_old
       lpd_prop <- llik_prop + lprior_prop
@@ -640,7 +633,7 @@ fit_vmmix <- function(data, ncomp, start_par = list(), method = "hmc", epsilon=0
         post.wt <- mem_p_univm(data.rad, par.mat.old, pi.mix.old, l.c.univm.old, ncores)
         clus.ind[ , iter] <- cID(post.wt, ncomp, runif(n.data))
         n.clus <- tabulate(clus.ind[ , iter], nbins = ncomp) #vector of component sizes
-        pi.mix.1 <- as.numeric(rdirichlet(1, (1 + n.clus))) #new mixture proportions
+        pi.mix.1 <- as.numeric(rdirichlet(1, (pmix.alpha + n.clus))) #new mixture proportions
         llik_new.pi <- llik_univm_full(data.rad, par.mat.old, pi.mix.1, l.c.univm.old, ncores)
       }
 
@@ -783,7 +776,7 @@ fit_vmmix <- function(data, ncomp, start_par = list(), method = "hmc", epsilon=0
         par.mat.prop <- q
         l.c.univm.prop <- log_const_univm_all(par.mat.prop)
 
-        lprior.prop <- sum(ldgamanum(q[1,], gam.loc, gam.scale))
+        lprior.prop <- sum((pmix.alpha-1) * log(pi.mix.1)) + sum(ldgamanum(q[1,], gam.loc, gam.scale))
 
         llik.prop <- llik_univm_full(data.rad, q, pi.mix.1, l.c.univm.prop, ncores)
 
@@ -854,7 +847,7 @@ fit_vmmix <- function(data, ncomp, start_par = list(), method = "hmc", epsilon=0
         post.wt <- mem_p_univm(data.rad, par.mat.old, pi.mix.old, l.c.univm.old, ncores)
         clus.ind[ , iter] <- cID(post.wt, ncomp, runif(n.data))
         n.clus <- tabulate(clus.ind[ , iter], nbins = ncomp) #vector of component sizes
-        pi.mix.1 <- as.numeric(rdirichlet(1, (1 + n.clus))) #new mixture proportions
+        pi.mix.1 <- as.numeric(rdirichlet(1, (pmix.alpha + n.clus))) #new mixture proportions
         llik_new.pi <- llik_univm_full(data.rad, par.mat.old, pi.mix.1, l.c.univm.old, ncores)
       }
 
@@ -874,7 +867,7 @@ fit_vmmix <- function(data, ncomp, start_par = list(), method = "hmc", epsilon=0
       lprior_old <- MC$lprior
 
       llik_prop <- llik_univm_full(data.rad, prop.mat, pi.mix.1, l.c.univm.prop, ncores)
-      lprior_prop <- sum(ldgamanum(k.1.prop, 0, 1000))
+      lprior_prop <- sum((pmix.alpha-1) * log(pi.mix.1)) + sum(ldgamanum(k.1.prop, gam.loc, gam.scale))
 
       lpd_old <- llik_old + lprior_old
       lpd_prop <- llik_prop + lprior_prop
@@ -978,6 +971,7 @@ fit_vmmix <- function(data, ncomp, start_par = list(), method = "hmc", epsilon=0
                  "epsilon.random" = epsilon.random, "epsilon" = epsilon_ave,
                  "L.random" = L.random, "L" = L_ave, "type" = "uni",
                  "propscale.final" = propscale_final, "data" = data.rad,
+                 "gam.loc" = gam.loc, "gam.scale" = gam.scale, "pmix.alpha" = pmix.alpha,
                  "n.data" = n.data, "ncomp" = ncomp, "n.iter" = n.iter)
 
   class(result) <- "angmcmc"