Export init and sample lm function for genMCMC example

Author: bking124

NAMESPACE

@@ -17,13 +17,15 @@ export(gbm_star)
 export(genEM_star)
 export(genMCMC_star)
 export(getEffSize)
+export(init_lm_gprior)
 export(lm_star)
 export(plot_coef)
 export(plot_fitted)
 export(plot_pmf)
 export(pvals)
 export(randomForest_star)
 export(round_floor)
+export(sample_lm_gprior)
 export(simBaS)
 export(simulate_nb_friedman)
 export(simulate_nb_lm)

R/STAR_Bayesian.R

@@ -280,15 +280,14 @@ blm_star <- function(y, X, X_test = NULL,
 #' are estimated using moments (means and variances) of \code{y}.
 #'
 #' @examples
-#' \donttest{
 #' # Simulate data with count-valued response y:
 #' sim_dat = simulate_nb_lm(n = 100, p = 5)
 #' y = sim_dat$y; X = sim_dat$X
 #'
 #' # STAR: log-transformation:
 #' fit_log = genMCMC_star(y = y,
-#'                          sample_params = function(y, params) rSTAR:::sample_lm_gprior(y, X, params),
-#'                          init_params = function(y) rSTAR:::init_lm_gprior(y, X),
+#'                          sample_params = function(y, params) sample_lm_gprior(y, X, params),
+#'                          init_params = function(y) init_lm_gprior(y, X),
 #'                          transformation = 'log')
 #' # Posterior mean of each coefficient:
 #' coef(fit_log)
@@ -297,14 +296,13 @@ blm_star <- function(y, X, X_test = NULL,
 #' fit_log$WAIC
 #'
 #' # MCMC diagnostics:
-#' plot(as.ts(fit_log$post.coefficients[,1:3]))
+#' plot(as.ts(fit_log$post.beta[,1:3]))
 #'
 #' # Posterior predictive check:
 #' hist(apply(fit_log$post.pred, 1,
 #'            function(x) mean(x==0)), main = 'Proportion of Zeros', xlab='');
 #' abline(v = mean(y==0), lwd=4, col ='blue')
 #'
-#'}
 #' @export
 genMCMC_star = function(y,
                      sample_params,
@@ -387,7 +385,7 @@ genMCMC_star = function(y,
     stop("The sample_params() function must return 'mu', 'sigma', and 'coefficients'")
 
   # Does the sampler return beta? If so, we want to store separately
-  beta_sampled = !is.null(params$coefficients$beta)
+  beta_sampled = !is.null(params$coefficients[["beta"]])
 
   #Does the sampler return mu_test
   testpoints = !is.null(params$mu_test)
@@ -533,7 +531,7 @@ genMCMC_star = function(y,
       if(nsi==1){
         print("Burn-In Period")
       } else if (nsi < nburn){
-        computeTimeRemaining(nsi, timer0, nstot, nrep = 4000)
+        computeTimeRemaining(nsi, timer0, nburn, nrep = 4000)
       } else if (nsi==nburn){
         print("Starting sampling")
         timer1 = proc.time()[3]
@@ -669,7 +667,7 @@ genMCMC_star = function(y,
 #' X_nonlin = as.matrix(X[,(1:3)])
 #'
 #' # STAR: nonparametric transformation
-#' fit <- bam_star(y,X_lin, X_nonlin)
+#' fit <- bam_star(y,X_lin, X_nonlin, nburn=1000, nskip=0)
 #'
 #' # Posterior mean of each coefficient:
 #' coef(fit)
@@ -819,8 +817,8 @@ bam_star = function(y, X_lin, X_nonlin, splinetype="orthogonal",
 #' y = sim_dat$y; X = sim_dat$X
 #'
 #' # BART-STAR with log-transformation:
-#' fit_log = bart_star(y = y, X = X,
-#'                          transformation = 'log', save_y_hat = TRUE)
+#' fit_log = bart_star(y = y, X = X, transformation = 'log',
+#'                     save_y_hat = TRUE, nburn=1000, nskip=0)
 #'
 #' # Fitted values
 #' plot_fitted(y = sim_dat$Ey,
@@ -1116,7 +1114,7 @@ bart_star = function(y,
       if(nsi==1){
         print("Burn-In Period")
       } else if (nsi < nburn){
-        computeTimeRemaining(nsi, timer0, nstot, nrep = 4000)
+        computeTimeRemaining(nsi, timer0, nburn, nrep = 4000)
       } else if (nsi==nburn){
         print("Starting sampling")
         timer1 = proc.time()[3]
@@ -1454,7 +1452,7 @@ spline_star = function(y,
       if(nsi==1){
         print("Burn-In Period")
       } else if (nsi < nburn){
-         computeTimeRemaining(nsi, timer0, nstot, nrep = 4000)
+         computeTimeRemaining(nsi, timer0, nburn, nrep = 4000)
       } else if (nsi==nburn){
         print("Starting sampling")
         timer1 = proc.time()[3]
@@ -1467,3 +1465,153 @@ spline_star = function(y,
 
   return(list(post_ytilde=post_ytilde))
 }
+
+#' Initialize linear regression parameters assuming a g-prior
+#'
+#' Initialize the parameters for a linear regression model assuming a
+#' g-prior for the coefficients.
+#'
+#' @param y \code{n x 1} vector of data
+#' @param X \code{n x p} matrix of predictors
+#' @param X_test \code{n0 x p} matrix of predictors at test points (default is NULL)
+#'
+#' @return a named list \code{params} containing at least
+#' \enumerate{
+#' \item \code{mu}: vector of conditional means (fitted values)
+#' \item \code{sigma}: the conditional standard deviation
+#' \item \code{coefficients}: a named list of parameters that determine \code{mu}
+#' }
+#' Additionally, if X_test is not NULL, then the list includes an element
+#' \code{mu_test}, the vector of conditional means at the test points
+#'
+#' @note The parameters in \code{coefficients} are:
+#' \itemize{
+#' \item \code{beta}: the \code{p x 1} vector of regression coefficients
+#' components of \code{beta}
+#' }
+#'
+#' @examples
+#' # Simulate data for illustration:
+#' sim_dat = simulate_nb_lm(n = 100, p = 5)
+#' y = sim_dat$y; X = sim_dat$X
+#'
+#' # Initialize:
+#' params = init_lm_gprior(y = y, X = X)
+#' names(params)
+#' names(params$coefficients)
+#'
+#' @export
+init_lm_gprior = function(y, X, X_test=NULL){
+
+  # Initialize the linear model:
+  n = nrow(X); p = ncol(X)
+
+  # Regression coefficients: depending on p >= n or p < n
+  if(p >= n){
+    beta = sampleFastGaussian(Phi = X, Ddiag = rep(1, p), alpha = y)
+  } else beta = lm(y ~ X - 1)$coef
+
+  # Fitted values:
+  mu = X%*%beta
+
+  #Mean at the test points (if passed in)
+  if(!is.null(X_test)) mu_test = X_test%*%beta
+
+  # Observation SD:
+  sigma = sd(y - mu)
+
+  # Named list of coefficients:
+  coefficients = list(beta = beta)
+
+  result = list(mu = mu, sigma = sigma, coefficients = coefficients)
+  if(!is.null(X_test)){
+    result = c(result, list(mu_test = mu_test))
+  }
+  return(result)
+}
+#' Sample the linear regression parameters assuming a g-prior
+#'
+#' Sample the parameters for a linear regression model assuming a
+#' g-prior for the  coefficients.
+#'
+#' @param y \code{n x 1} vector of data
+#' @param X \code{n x p} matrix of predictors
+#' @param params the named list of parameters containing
+#' \enumerate{
+#' \item \code{mu}: vector of conditional means (fitted values)
+#' \item \code{sigma}: the conditional standard deviation
+#' \item \code{coefficients}: a named list of parameters that determine \code{mu}
+#' }
+#' @param psi the prior variance for the g-prior
+#' @param XtX the \code{p x p} matrix of \code{crossprod(X)} (one-time cost);
+#' if NULL, compute within the function
+#' @param X_test matrix of predictors at test points (default is NULL)
+#'
+#' @return The updated named list \code{params} with draws from the full conditional distributions
+#' of \code{sigma} and \code{coefficients} (along with updated \code{mu} and \code{mu_test} if applicable).
+#'
+#' @note The parameters in \code{coefficients} are:
+#' \itemize{
+#' \item \code{beta}: the \code{p x 1} vector of regression coefficients
+#' components of \code{beta}
+#' }
+#'
+#' @examples
+#' # Simulate data for illustration:
+#' sim_dat = simulate_nb_lm(n = 100, p = 5)
+#' y = sim_dat$y; X = sim_dat$X
+#' # Initialize:
+#' params = init_lm_gprior(y = y, X = X)
+#' # Sample:
+#' params = sample_lm_gprior(y = y, X = X, params = params)
+#' names(params)
+#' names(params$coefficients)
+#'
+#' @import truncdist
+#' @export
+sample_lm_gprior = function(y, X, params, psi = NULL, XtX = NULL, X_test=NULL){
+
+  # Dimensions:
+  n = nrow(X); p = ncol(X)
+
+  if(is.null(psi)) psi = n # default
+
+  # For faster computations:
+  if(is.null(XtX)) XtX = crossprod(X)
+
+  # Access elements of the named list:
+  sigma = params$sigma  # Observation SD
+  coefficients = params$coefficients # Coefficients to access below:
+
+  beta = coefficients$beta;              # Regression coefficients (including intercept)
+
+  # Sample the regression coefficients:
+  Q_beta = 1/sigma^2*(1+psi)/(psi)*XtX
+  ell_beta = 1/sigma^2*crossprod(X, y)
+  ch_Q = chol(Q_beta)
+  beta = backsolve(ch_Q,
+                   forwardsolve(t(ch_Q), ell_beta) +
+                     rnorm(p))
+
+  # Conditional mean:
+  mu = X%*%beta
+
+  #Mean at the test points (if passed in)
+  if(!is.null(X_test)) mu_test = X_test%*%beta
+
+  # Observation SD:
+  sigma =  1/sqrt(rgamma(n = 1,
+                         shape = .001 + n/2,
+                         rate = .001 + sum((y - mu)^2)/2))
+
+  # Update the coefficients:
+  coefficients$beta = beta
+
+  result = list(mu = mu, sigma = sigma, coefficients = coefficients)
+  if(!is.null(X_test)){
+    result = c(result, list(mu_test = mu_test))
+  }
+  return(result)
+}
+
+

_pkgdown.yml

@@ -26,6 +26,8 @@ reference:
   - bart_star
   - spline_star
   - genMCMC_star
+  - init_lm_gprior
+  - sample_lm_gprior
 - title: WarpDLMs (Count Time Series Modeling)
   desc: Methods for warped Dynamic Linear Models to perform Bayesian inference and forecasting for count-valued time series
 - contents: