updated simulations, plot and confint fns

R/bw_table.R

@@ -6,13 +6,14 @@ num_sims = 100
 y = 0.0
 x = 0.0
 #true density function
-true_dens = stats::dnorm(y, mean = 0, sd=1)
+true_dens = stats::dnorm(y, mean = x, sd=1)
+
+hmse=0.48
 
 #grid for bandwidth range
-h_grid = c(0.5, 0.6, 0.7, 0.8, 0.9, 1, 1.1, 1.2, 1.3, 1.4, 1.5)*0.5
+h_grid = c(0.5, 0.6, 0.7, 0.8, 0.9, 1, 1.1, 1.2, 1.3, 1.4, 1.5)*hmse
 
 bw_testing = function(s){
-  print(s)
   # simulating data
   x_data = matrix(rnorm(n, mean=0, sd=1))
   y_data = matrix(rnorm(n, mean=x_data, sd=1))
@@ -44,9 +45,11 @@ bw_testing = function(s){
 
 table_data = matrix(ncol=8)
 # running simulations
+RNGkind(kind = "L'Ecuyer-CMRG")
+set.seed(42)
 for (i in 1:length(h_grid)){
  h = h_grid[i]
- object = parallel::mclapply(1:num_sims, bw_testing, mc.cores =8)
+ object = parallel::mclapply(1:num_sims, bw_testing, mc.cores=8)
   int_avg = matrix(0L, nrow = num_sims, ncol = 8)
   for (l in 1:num_sims){
     int_avg[l,] = object[[l]][1,]
@@ -57,6 +60,6 @@ table_data = table_data[2:nrow(table_data), ]
 table_data[, 2] = abs(table_data[,2])
 colnames(table_data) = c("BW", "bias", "sd", "rmse", "WBC CE", "RBC CE",
                          "WBC AL", "RBC AL")
-rownames(table_data) = h_grid/0.5
+rownames(table_data) = h_grid/hmse
 #print table
 table_data

R/lpcde/R/lpcde_methods.R

@@ -321,10 +321,8 @@ vcov.lpcde = function(object, ...) {
 #'
 #' @param object Class "lpdensity" object, obtained by calling \code{\link{lpcde}}.
 #' @param parm Integer, indicating which parameters are to be given confidence intervals.
-#' @param level Numeric scalar between 0 and 1, the significance level for computing
-#' confidence intervals
-#' @param alpha Numeric scalar between 0 and 1, specifies the significance level for plotting
-#'   confidence intervals/bands.
+#' @param level Numeric scalar between 0 and 1, the confidence level for computing
+#' confidence intervals/bands. Equivalent to (1-significance level).
 #' @param CIuniform \code{TRUE} or \code{FALSE} (default), plotting either pointwise confidence intervals (\code{FALSE}) or
 #'   uniform confidence bands (\code{TRUE}).
 #' @param CIsimul Positive integer, specifies the number of simulations used to construct critical values (default is \code{2000}). This
@@ -368,17 +366,15 @@ vcov.lpcde = function(object, ...) {
 #' confint(model1)
 #'
 #' @export
-confint.lpcde <- function(object, parm = NULL, level = NULL, CIuniform=FALSE, CIsimul=2000, alpha=0.05, ...){
+confint.lpcde <- function(object, parm = NULL, level=0.95, CIuniform=FALSE, CIsimul=2000, ...){
   x = object
   if (class(x)[1] == "lpbwcde") {
     stop("The confint method does not support \"lpbwcde\" objects.\n")
   }
   args = list(...)
+  alpha = 1-level
 
   if (!is.null(parm)) { args[['grid']] = parm }
-  if (!is.null(level)) { args[['alpha']] = 1 - level }
-
-  if (is.null(args[['alpha']])) { alpha = 0.05 } else { alpha = args[['alpha']] }
   if (is.null(args[['CIuniform']]))   { CIuniform = FALSE } else { CIuniform = args[['CIuniform']] }
   if (is.null(args[['CIsimul']]))   { CIsimul = 2000 } else { sep = args[['CIsimul']] }
 
@@ -786,8 +782,6 @@ plot.lpcde = function(..., alpha=NULL,type=NULL, lty=NULL, lwd=NULL, lcol=NULL,
     if (legend_default) {
       data_x$Sname = paste("Series", i, sep=" ")
       legendGroups = c(legendGroups, data_x$Sname)
-    } else {
-      data_x$Sname = legendGroups[i]
     }
 
     ########################################
@@ -875,7 +869,7 @@ plot.lpcde = function(..., alpha=NULL,type=NULL, lty=NULL, lwd=NULL, lcol=NULL,
   }
 
   if (is.null(xlabel)) {
-    xlabel <- ""
+    xlabel <- "y"
   }
 
   if (is.null(title)) {
@@ -885,6 +879,10 @@ plot.lpcde = function(..., alpha=NULL,type=NULL, lty=NULL, lwd=NULL, lcol=NULL,
   temp_plot <- temp_plot + ggplot2::labs(x=xlabel, y=ylabel) + ggplot2::ggtitle(title) +
     ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5))
 
+  if (nfig ==1){
+    temp_plot <- temp_plot + ggplot2::theme(legend.position = "none")
+  }
+
   # check plotting range vs estimation range
   if (!is.null(y_grid)) {
     if (min(y_grid) < estRangeL | max(y_grid) > estRangeR) {

R/lpcde_comparison.R

@@ -12,7 +12,7 @@ y = as.matrix(iris[, 1], ncol=1)
 x = as.matrix(iris[, 3], ncol=1)
 x_evals = quantile(x)
 
-#scatter plot of data with chosen conditioning values
+#scatter plot of data with chosen conditioning values (See Fig 3)
 plot(x,y, xlab = "Petal length", ylab = "Sepal length")
 abline(v=x_evals[2], col="green")
 abline(v=x_evals[3], col="green")
@@ -37,27 +37,33 @@ npbws_q3 = npcdensbw(x, y)
 np_est_q3 = fitted(npcdens(exdat=rep(x_evals[4], 30), eydat=lpcde_est_q3$Estimate[,1], bws=npbws_q3))
 
 #plotting all CDEs
+# Figure 4(a)
 plot(1, main=TeX(r'($f(y|x=1.6)$)'),  xlab="Sepal length", ylab="density", ylim=c(-0.8, 1.6), xlim=c(4.8, 6.9), cex=1.5)
 lines(lpcde_est_q1$Estimate[, 1], lpcde_est_q1$Estimate[,3])
 lines(hdr_est_q1$y, hdr_est_q1$z, col=2)
 lines(lpcde_est_q1$Estimate[, 1], np_est_q1, col=3)
 legend('topright',lwd=1, legend=c('lpcde', 'hdrcde', 'np'), col=c(1,2,3))
 
+# Figure 5(a)
 plot(1, main=TeX(r'($f(y|x=4.35)$)'),  xlab="Sepal length", ylab="density", ylim=c(-0.8, 1.6), xlim=c(4.8, 6.9))
 lines(lpcde_est_q2$Estimate[, 1], lpcde_est_q2$Estimate[,3])
 lines(hdr_est_q2$y, hdr_est_q2$z, col=2)
 lines(lpcde_est_q2$Estimate[, 1], np_est_q2, col=3)
 legend('topright',lwd=1, legend=c('lpcde', 'hdrcde', 'np'), col=c(1,2,3))
 
+# Figure 6(a)
 plot(1, main=TeX(r'($f(y|x=5.1)$)'),  xlab="Sepal length", ylab="density", ylim=c(-0.8, 1.6), xlim=c(4.8, 6.9))
 lines(lpcde_est_q3$Estimate[, 1], lpcde_est_q3$Estimate[,3])
 lines(hdr_est_q3$y, hdr_est_q3$z, col=2)
 lines(lpcde_est_q3$Estimate[, 1], np_est_q3, col=3)
 legend('topright',lwd=1, legend=c('lpcde', 'hdrcde', 'np'), col=c(1,2,3))
 
-# lpcde estimates with CI
-plot(lpcde_est_q1, xlabel="Sepal length", ylabel="density", title=TeX(r'($f(y|x=1.6)$ with confidence bands)')) + ggplot2::ylim(-0.8, 1.6) + ggplot2::xlim(4.8, 6.9) + ggplot2::theme(legend.position='none', panel.grid.major = ggplot2::element_blank(), panel.grid.minor = ggplot2::element_blank(), panel.background = ggplot2::element_blank(), text = ggplot2::element_text(size = 7))
+# Figure 4(b)
+plot(lpcde_est_q1, xlabel="Sepal length", ylabel="density", title=TeX(r'($f(y|x=1.6)$ with confidence bands)'))+ ggplot2::ylim(-0.8, 1.6) + ggplot2::xlim(4.8, 6.9) + ggplot2::theme(legend.position='none', panel.grid.major = ggplot2::element_blank(), panel.grid.minor = ggplot2::element_blank(),
+                                                                                                                                                                                     panel.background = ggplot2::element_blank(), text = ggplot2::element_text(size = 7))
+# Figure 5(b)
 plot(lpcde_est_q2,xlabel="Sepal length", ylabel="density", title=TeX(r'($f(y|x=4.35)$ with confidence bands)'))+ ggplot2::ylim(-0.8, 1.6) + ggplot2::xlim(4.8, 6.9) + ggplot2::theme(legend.position='none', panel.grid.major = ggplot2::element_blank(), panel.grid.minor = ggplot2::element_blank(),
                                                                                                                                                                                      panel.background = ggplot2::element_blank(), text = ggplot2::element_text(size = 7))
+# Figure 6(b)
 plot(lpcde_est_q3, xlabel="Sepal length", ylabel="density", title=TeX(r'($f(y|x=5.1)$ with confidence bands)'))+ ggplot2::ylim(-0.8, 1.6) + ggplot2::xlim(4.8, 6.9) + ggplot2::theme(legend.position='none', panel.grid.major = ggplot2::element_blank(), panel.grid.minor = ggplot2::element_blank(),
                                                                                                                                                                                      panel.background = ggplot2::element_blank(), text = ggplot2::element_text(size = 7))

R/lpcde_replication.R

@@ -4,7 +4,7 @@
 library("lpcde")
 
 #Setting up simulation
-set.seed(42)
+set.seed(30)
 n=1000
 x_data = matrix(rnorm(n, mean=0, sd=1))
 y_data = matrix(rnorm(n, mean=x_data, sd=1))

R/pw_table.R

@@ -6,8 +6,9 @@ num_sims = 100
 y_points = c(0.0, 0.8, 1.5)
 x_points = c(0.0, 0.8, 1.5)
 
+set.seed(42)
+RNGkind(kind = "L'Ecuyer-CMRG")
 pw_testing = function(s){
-  print(s)
   # simulating data
   x_data = matrix(rnorm(n, mean=0, sd=1))
   y_data = matrix(rnorm(n, mean=x_data, sd=1))
@@ -39,7 +40,7 @@ pw_testing = function(s){
 }
 
 pw_table = list()
-bw_mat = matrix(c(0.51, 0.58, 0.8, 0.6, 0.6, 0.6, 1, 0.8, .9), nrow=3)
+bw_mat = matrix(c(0.5, 0.58, 0.82, 0.6, 0.6, 0.6, 1, 0.8, 0.8), nrow=3)
 # running simulations
 for (i in 1:length(x_points)){
   x = x_points[i]