Merge branch 'main' into develop

.github/workflows/CI.yml

@@ -29,3 +29,4 @@ jobs:
           files: ./R/lpcde/cobertura.xml
           fail_ci_if_error: true
           verbose: true
+          version: "v0.1.15"

.github/workflows/draft-pdf.yml

@@ -0,0 +1,29 @@
+name: Draft PDF
+on:
+  push:
+    paths:
+      - paper.md
+      - paper.bib
+      - .github/workflows/draft-pdf.yml
+
+jobs:
+  paper:
+    runs-on: ubuntu-latest
+    name: Paper Draft
+    steps:
+      - name: Checkout
+        uses: actions/checkout@v4
+      - name: Build draft PDF
+        uses: openjournals/openjournals-draft-action@master
+        with:
+          journal: joss
+          # This should be the path to the paper within your repo.
+          paper-path: paper.md
+      - name: Upload
+        uses: actions/upload-artifact@v4
+        with:
+          name: paper
+          # This is the output path where Pandoc will write the compiled
+          # PDF. Note, this should be the same directory as the input
+          # paper.md
+          path: paper.pdf

.gitignore

@@ -7,6 +7,13 @@ R/lpcde/.Rhistory
 #dont include tar
 R/*.tar.gz
 
+#replication files
+*.tex
+*.txt
+
 *.html
 .Rproj.user
 */.DS_Store
+
+#emacs
+auto/

R/bw_table.R

@@ -1,18 +1,19 @@
 # Script for generating table 1 in Cattaneo, Chandak, Jansson and Ma (2022)
 
 #parameters
-n = 1000
-num_sims = 1000
+n = 2000
+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.6
+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,20 +45,24 @@ 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)
+  h = h_grid[i]
+  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,]
   }
-table_data = rbind(table_data, (colSums(int_avg)/num_sims))
+  table_data = rbind(table_data, (colSums(int_avg)/num_sims))
 }
 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.6
+rownames(table_data) = h_grid/hmse
 #print table
+library(xtable)
 table_data
-#print(xtable(table_data), file='bw_table.tex')
+print(xtable(table_data), file="bw_table.tex")
+

R/lpcde/DESCRIPTION

@@ -1,13 +1,13 @@
 Package: lpcde
 Type: Package
 Title: Boundary Adaptive Local Polynomial Conditional Density Estimator
-Version: 0.1.2
+Version: 0.1.4
 Authors@R: person(given = "Rajita",
            family = "Chandak",
            role = c("aut", "cre"),
            email = "rchandak@princeton.edu" )
 Maintainer: Rajita Chandak <rchandak@princeton.edu>
-Description: Tools for estimation and inference of conditional densities, derivatives and functions. This is the companion software for Cattaneo, Chandak, Jansson and Ma (2024) <arxiv:2204.10359>.
+Description: Tools for estimation and inference of conditional densities, derivatives and functions. This is the companion software for Cattaneo, Chandak, Jansson and Ma (2024) <doi:10.48550/arXiv.2204.10359>.
 Depends: R (>= 3.3.0)
 License: GPL-2
 Encoding: UTF-8

R/lpcde/NEWS.md

@@ -1,3 +1,12 @@
+# lpcde 0.1.4
+- Cleaned up standard plot output
+- Consolidated level and alpha parameters for confint
+
+# lpcde 0.1.3
+- Fixed uniform CB construction
+- Added rbc flag for plotting
+- Fixed re-scaling of inputs
+
 # lpcde 0.1.2
 - Added flags for returning normalized and/or non-negative estimate
 - Added option to change automatic grid spacing (options include quantile and equally-spaced grid generation). Functionality available for lpcde and lpbwcde.

R/lpcde/R/RcppExports.R

@@ -4,7 +4,7 @@
 #' @title All Sums in C++ (Internal Function)
 #' @description Function that prints all combinations of natural numbers that
 #' add up to target value.
-#' @param Target target value for sum.
+#' @param target Target value for sum.
 #' @return List of combinations that add up to target value.
 #' @keywords internal
 print_all_sumC <- function(target) {

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/src/utilsC.cpp

@@ -52,7 +52,7 @@ void print_all_sum_rec(
 //' @title All Sums in C++ (Internal Function)
 //' @description Function that prints all combinations of natural numbers that
 //' add up to target value.
-//' @param Target target value for sum.
+//' @param target Target value for sum.
 //' @return List of combinations that add up to target value.
 //' @keywords internal
 // [[Rcpp::export]]

R/lpcde_comparison.R

@@ -12,8 +12,8 @@ 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
-plot(x,y)
+#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")
 abline(v=x_evals[4], col="green")
@@ -37,48 +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
-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))
-lines(lpcde_est_q1$Estimate[, 1], lpcde_est_q1$Estimate[,3])
+# 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')
-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')
-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')
-
-
-
-
-
-#open3d()
-#
-#lines3d(x_evals[2], lpcde_est_q1$Estimate[,1], lpcde_est_q1$Estimate[,3])
-#lines3d(x_evals[3], lpcde_est_q2$Estimate[,1], lpcde_est_q2$Estimate[,3])
-#lines3d(x_evals[4], lpcde_est_q3$Estimate[,1], lpcde_est_q3$Estimate[,3])
-#
-#lines3d(x_evals[2], hdr_est_q1$y, hdr_est_q1$z, color='red')
-#lines3d(x_evals[3], hdr_est_q2$y, hdr_est_q2$z, color='red')
-#lines3d(x_evals[4], hdr_est_q3$y, hdr_est_q3$z, color='red')
-#
-#lines3d(x_evals[2], lpcde_est_q1$Estimate[,1], np_est_q1, color='blue')
-#lines3d(x_evals[3], lpcde_est_q2$Estimate[,1], np_est_q2, color='blue')
-#lines3d(x_evals[4], lpcde_est_q3$Estimate[,1], np_est_q3, color='blue')
-#
-#axes3d()
-#title3d('', '', 'Petal Length', 'Sepal Length', 'density')
-#legend3d("topright", c("lpcde", "hdrcde", "np"), pch = c(1,25, 16), col = c('black', 'red', 'blue'), cex=1)
-
+# 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))
@@ -31,5 +31,5 @@ lines(model1$Estimate[,1], dnorm(y_grid), col=3)
 legend('topleft',lwd=1, legend=c('standard estimate', 'regularized estimate', 'true density'), col=c(1,2,3))
 
 #A simple plot (see Fig 2)
-plot(model1, CIuniform = TRUE, rbc=TRUE, xlabel="y") + ggplot2::theme(legend.position = "none")
+plot(model1, CIuniform = TRUE, rbc=TRUE, xlabel="y")