Automated build.

.setup/cpp/cTMed-mc-phi-sigma.cpp

@@ -0,0 +1,60 @@
+// -----------------------------------------------------------------------------
+// edit .setup/cpp/cTMed-mc-phi-sigma.cpp
+// Ivan Jacob Agaloos Pesigan
+// -----------------------------------------------------------------------------
+
+#include <RcppArmadillo.h>
+// [[Rcpp::depends(RcppArmadillo)]]
+// [[Rcpp::export(.MCPhiSigma)]]
+Rcpp::List MCPhiSigma(const arma::vec& theta, const arma::mat& vcov_theta,
+                      const arma::uword& R, bool test_phi = true) {
+  Rcpp::List output(R);
+  arma::uword n = theta.n_elem;
+  arma::uword p = (-1 + std::sqrt(1 + 24 * n)) / 6;
+  arma::uword q = (p * (p + 1)) / 2;
+  arma::uword index = 0;
+  arma::vec v_i(n, arma::fill::none);
+  arma::mat phi_i(p, p, arma::fill::none);
+  arma::vec phi_vec_i(p * p, arma::fill::none);
+  arma::mat sigma_i(p, p, arma::fill::none);
+  arma::vec sigma_vech_i(q, arma::fill::none);
+  arma::vec eigval;
+  arma::mat eigvec;
+  for (arma::uword i = 0; i < R; i++) {
+    bool run = true;
+    Rcpp::List output_i(2);
+    while (run) {
+      // generate data
+      v_i = arma::mvnrnd(theta, vcov_theta);
+      phi_vec_i = v_i(arma::span(0, (p * p) - 1));
+      sigma_vech_i = v_i(arma::span(p * p, n - 1));
+      // test phi
+      phi_i = arma::reshape(phi_vec_i, p, p);
+      if (test_phi) {
+        if (TestPhi(phi_i)) {
+          run = false;
+        }
+      } else {
+        run = false;
+      }
+      if (run == false) {
+        // test sigma
+        index = 0;
+        for (arma::uword i = 0; i < p; ++i) {
+          for (arma::uword j = i; j < p; ++j) {
+            sigma_i(i, j) = sigma_vech_i(index);
+            sigma_i(j, i) = sigma_vech_i(index);
+            index++;
+          }
+        }
+        arma::eig_sym(eigval, eigvec, sigma_i);
+        eigval.transform([](double val) { return std::max(val, 1e-8); });
+        sigma_i = eigvec * arma::diagmat(eigval) * eigvec.t();
+      }
+    }
+    output_i[0] = phi_i;
+    output_i[1] = sigma_i;
+    output[i] = output_i;
+  }
+  return output;
+}

.setup/latex/pdf/.gitignore

@@ -0,0 +1,5 @@
+*
+*/
+!*.pdf
+!.gitignore
+!bib.bib

NAMESPACE

@@ -44,6 +44,7 @@ export(MCIndirectCentral)
 export(MCMed)
 export(MCMedStd)
 export(MCPhi)
+export(MCPhiSigma)
 export(MCTotalCentral)
 export(Med)
 export(MedStd)

NEWS.md

@@ -3,6 +3,8 @@
 ## Patch
 
 * Latest development version.
+* Added the `MCPhiSigma()` function.
+* Minor edits to tests.
 
 # cTMed 1.0.4
 

R/RcppExports.R

@@ -45,6 +45,10 @@
     .Call(`_cTMed_MCPhiSigmaI`, theta, vcov_theta, test_phi)
 }
 
+.MCPhiSigma <- function(theta, vcov_theta, R, test_phi = TRUE) {
+    .Call(`_cTMed_MCPhiSigma`, theta, vcov_theta, R, test_phi)
+}
+
 .MCPhi <- function(phi, vcov_phi_vec_l, R, test_phi = TRUE) {
     .Call(`_cTMed_MCPhi`, phi, vcov_phi_vec_l, R, test_phi)
 }

R/cTMed-mc-med-std.R

@@ -80,12 +80,9 @@
 #'
 #' @author Ivan Jacob Agaloos Pesigan
 #'
-#' @param vcov_theta Numeric matrix.
-#'   The sampling variance-covariance matrix of
-#'   \eqn{\mathrm{vec} \left( \boldsymbol{\Phi} \right)} and
-#'   \eqn{\mathrm{vech} \left( \boldsymbol{\Sigma} \right)}
 #' @inheritParams IndirectStd
 #' @inheritParams MCMed
+#' @inheritParams MCPhiSigma
 #' @inherit IndirectStd references
 #'
 #' @return Returns an object

R/cTMed-mc-phi-sigma.R

@@ -0,0 +1,214 @@
+#' Generate Random Drift Matrices
+#' and Process Noise Covariance Matrices
+#' Using the Monte Carlo Method
+#'
+#' This function generates random
+#' drift matrices \eqn{\boldsymbol{\Phi}}
+#' and process noise covariabces matrices \eqn{\boldsymbol{\Sigma}}
+#' using the Monte Carlo method.
+#'
+#' @details
+#' ## Monte Carlo Method
+#'   Let \eqn{\boldsymbol{\theta}} be
+#'   a vector that combines
+#'   \eqn{\mathrm{vec} \left( \boldsymbol{\Phi} \right)},
+#'   that is,
+#'   the elements of the \eqn{\boldsymbol{\Phi}} matrix
+#'   in vector form sorted column-wise and
+#'   \eqn{\mathrm{vech} \left( \boldsymbol{\Sigma} \right)},
+#'   that is,
+#'   the unique elements of the \eqn{\boldsymbol{\Sigma}} matrix
+#'   in vector form sorted column-wise.
+#'   Let \eqn{\hat{\boldsymbol{\theta}}} be
+#'   a vector that combines
+#'   \eqn{\mathrm{vec} \left( \hat{\boldsymbol{\Phi}} \right)} and
+#'   \eqn{\mathrm{vech} \left( \hat{\boldsymbol{\Sigma}} \right)}.
+#'   Based on the asymptotic properties of maximum likelihood estimators,
+#'   we can assume that estimators are normally distributed
+#'   around the population parameters.
+#'   \deqn{
+#'   	\hat{\boldsymbol{\theta}}
+#'   	\sim
+#'   	\mathcal{N}
+#'   	\left(
+#'   	\boldsymbol{\theta},
+#'   	\mathbb{V} \left( \hat{\boldsymbol{\theta}} \right)
+#'   	\right)
+#'   }
+#'   Using this distributional assumption,
+#'   a sampling distribution of \eqn{\hat{\boldsymbol{\theta}}}
+#'   which we refer to as \eqn{\hat{\boldsymbol{\theta}}^{\ast}}
+#'   can be generated by replacing the population parameters
+#'   with sample estimates,
+#'   that is,
+#'   \deqn{
+#'   	\hat{\boldsymbol{\theta}}^{\ast}
+#'   	\sim
+#'   	\mathcal{N}
+#'   	\left(
+#'   	\hat{\boldsymbol{\theta}},
+#'   	\hat{\mathbb{V}} \left( \hat{\boldsymbol{\theta}} \right)
+#'   	\right) .
+#'   }
+#'
+#' @author Ivan Jacob Agaloos Pesigan
+#'
+#' @inheritParams IndirectStd
+#' @inheritParams MCPhi
+#' @param vcov_theta Numeric matrix.
+#'   The sampling variance-covariance matrix of
+#'   \eqn{\mathrm{vec} \left( \boldsymbol{\Phi} \right)} and
+#'   \eqn{\mathrm{vech} \left( \boldsymbol{\Sigma} \right)}
+#'
+#' @return Returns an object
+#'   of class `ctmedmc` which is a list with the following elements:
+#'   \describe{
+#'     \item{call}{Function call.}
+#'     \item{args}{Function arguments.}
+#'     \item{fun}{Function used ("MCPhiSigma").}
+#'     \item{output}{A list simulated drift matrices.}
+#'   }
+#'
+#' @examples
+#' set.seed(42)
+#' phi <- matrix(
+#'   data = c(
+#'     -0.357, 0.771, -0.450,
+#'     0.0, -0.511, 0.729,
+#'     0, 0, -0.693
+#'   ),
+#'   nrow = 3
+#' )
+#' colnames(phi) <- rownames(phi) <- c("x", "m", "y")
+#' sigma <- matrix(
+#'   data = c(
+#'     0.24455556, 0.02201587, -0.05004762,
+#'     0.02201587, 0.07067800, 0.01539456,
+#'     -0.05004762, 0.01539456, 0.07553061
+#'   ),
+#'   nrow = 3
+#' )
+#' MCPhiSigma(
+#'   phi = phi,
+#'   sigma = sigma,
+#'   vcov_theta = 0.1 * diag(15),
+#'   R = 100L # use a large value for R in actual research
+#' )
+#'
+#' @family Continuous Time Mediation Functions
+#' @keywords cTMed mc
+#' @export
+MCPhiSigma <- function(phi,
+                       sigma,
+                       vcov_theta,
+                       R,
+                       test_phi = TRUE,
+                       ncores = NULL,
+                       seed = NULL) {
+  idx <- rownames(phi)
+  stopifnot(
+    idx == colnames(phi)
+  )
+  args <- list(
+    phi = phi,
+    sigma = sigma,
+    vcov_theta = vcov_theta,
+    R = R,
+    test_phi = test_phi,
+    ncores = ncores,
+    seed = seed,
+    method = "mc"
+  )
+  # nocov start
+  par <- FALSE
+  if (!is.null(ncores)) {
+    ncores <- as.integer(ncores)
+    if (ncores > R) {
+      ncores <- R
+    }
+    if (ncores > 1) {
+      par <- TRUE
+    }
+  }
+  if (par) {
+    os_type <- Sys.info()["sysname"]
+    if (os_type == "Darwin") {
+      fork <- TRUE
+    } else if (os_type == "Linux") {
+      fork <- TRUE
+    } else {
+      fork <- FALSE
+    }
+    if (fork) {
+      if (!is.null(seed)) {
+        set.seed(seed)
+      }
+      output <- parallel::mclapply(
+        X = seq_len(R),
+        FUN = function(i) {
+          return(
+            .MCPhiSigmaI(
+              theta = c(.Vec(phi), .Vech(sigma)),
+              vcov_theta = vcov_theta,
+              test_phi = test_phi
+            )
+          )
+        },
+        mc.cores = ncores
+      )
+    } else {
+      cl <- parallel::makeCluster(ncores)
+      on.exit(
+        parallel::stopCluster(cl = cl)
+      )
+      if (!is.null(seed)) {
+        parallel::clusterSetRNGStream(
+          cl = cl,
+          iseed = seed
+        )
+      }
+      output <- parallel::parLapply(
+        cl = cl,
+        X = seq_len(R),
+        fun = function(i) {
+          return(
+            .MCPhiSigmaI(
+              theta = c(.Vec(phi), .Vech(sigma)),
+              vcov_theta = vcov_theta,
+              test_phi = test_phi
+            )
+          )
+        }
+      )
+    }
+    # nocov end
+  } else {
+    if (!is.null(seed)) {
+      set.seed(seed)
+    }
+    output <- .MCPhiSigma(
+      theta = c(.Vec(phi), .Vech(sigma)),
+      vcov_theta = vcov_theta,
+      R = R,
+      test_phi = test_phi
+    )
+  }
+  output <- lapply(
+    X = output,
+    FUN = function(x) {
+      colnames(x[[1]]) <- rownames(x[[1]]) <- idx
+      return(x)
+    }
+  )
+  out <- list(
+    call = match.call(),
+    args = args,
+    fun = "MCPhiSigma",
+    output = output
+  )
+  class(out) <- c(
+    "ctmedmcphi",
+    class(out)
+  )
+  return(out)
+}

R/cTMed-methods-ctmedmcphi.R

@@ -30,11 +30,31 @@
 print.ctmedmcphi <- function(x,
                              digits = 4,
                              ...) {
-  base::print(
-    lapply(
-      X = x$output,
-      FUN = round,
-      digits = digits
+  if (x$fun == "MCPhi") {
+    base::print(
+      lapply(
+        X = x$output,
+        FUN = round,
+        digits = digits
+      )
     )
-  )
+  }
+  if (x$fun == "MCPhiSigma") {
+    base::print(
+      lapply(
+        X = x$output,
+        FUN = function(x) {
+          x[[1]] <- round(
+            x[[1]],
+            digits = digits
+          )
+          x[[2]] <- round(
+            x[[2]],
+            digits = digits
+          )
+          return(x)
+        }
+      )
+    )
+  }
 }

README.md

@@ -1,7 +1,7 @@
 cTMed
 ================
 Ivan Jacob Agaloos Pesigan
-2025-01-09
+2025-01-11
 
 <!-- README.md is generated from README.Rmd. Please edit that file -->