fixed bw error

nppackages · Oct 14, 2024 · 90c55d3 · 90c55d3
1 parent a642f48
commit 90c55d3
Show file tree

Hide file tree

Showing 7 changed files with 32 additions and 24 deletions.
diff --git a/R/lpcde/R/lpbwcde.R b/R/lpcde/R/lpbwcde.R
@@ -73,7 +73,7 @@
 #'
 lpbwcde <- function(y_data, x_data, x, y_grid=NULL, p=NULL, q=NULL, grid_spacing="", ng=NULL,
                     mu=NULL, nu=NULL, kernel_type=c("epanechnikov", "triangular", "uniform"),
-                    bw_type=c("mse-rot", "imse-rot"), regularize=NULL){
+                    bw_type=c("imse-rot", "mse-rot"), regularize=NULL){
   ################################################################################
   # Error Checking
   ################################################################################
@@ -262,7 +262,7 @@ lpbwcde <- function(y_data, x_data, x, y_grid=NULL, p=NULL, q=NULL, grid_spacing
   }
 
   Result = list(BW=BW,
-                 opt=list(x=x, p=p, q=q, mu=mu, nu=nu, kernel_type=kernel_type, n=n, ng=ng,
+                 opt=list(x=(x*sd_x + mx), p=p, q=q, mu=mu, nu=nu, kernel_type=kernel_type, n=n, ng=ng,
                           bw_type=bw_type,
                           data_min=min(y_data), data_max=max(y_data),
                           grid_min=min(y_grid), grid_max=max(y_grid)))

diff --git a/R/lpcde/R/lpbwcde_fns.R b/R/lpcde/R/lpbwcde_fns.R
@@ -90,8 +90,7 @@ bw_rot = function(y_data, x_data, y_grid, x, p, q, mu, nu, kernel_type, regulari
     } else{
       alpha = d + 2*min(p, q) + 2*max(mu, nu) + 1
     }
-    h = (v_dgp/bias_dgp[, 3])^(1/alpha)*n^(-1/alpha)
-    h = sd_y*sd_x*h
+    h = (abs(v_dgp/bias_dgp[, 3]))^(1/alpha)*n^(-1/alpha)
 
     if (regularize == TRUE){
       for (j in 1:ng){
@@ -275,7 +274,7 @@ bw_irot = function(y_data, x_data, y_grid, x, p, q, mu, nu, kernel_type, regular
     } else{
       alpha = d + 2*min(p, q) + 2*max(mu, nu) + 1
     }
-    h = (mean(v_dgp)/(2*mean(bias_dgp[, 3])))^(1/alpha)*n^(-1/alpha)
+    h = (abs(mean(v_dgp)/(2*mean(bias_dgp[, 3]))))^(1/alpha)*n^(-1/alpha)
     h = sd_y*sd_x*h
 
     if (regularize == TRUE){

diff --git a/R/lpcde/R/lpcde.R b/R/lpcde/R/lpcde.R
@@ -258,7 +258,7 @@ lpcde = function(x_data, y_data, y_grid=NULL, x=NULL, bw=NULL, p=NULL, q=NULL,
   # bw
   if (length(bw) == 0) {
     if (length(bw_type) == 0) {
-      bw_type = "mse-rot"
+      bw_type = "imse-rot"
       bw = lpbwcde(y_data=y_data, x_data=x_data, x=x, y_grid=y_grid, p=p, q=q, mu=mu,
                    nu=nu, kernel_type=kernel_type, bw_type=bw_type)$BW[,2]
     } else {

diff --git a/R/lpcde/R/lpcde_fns.R b/R/lpcde/R/lpcde_fns.R
@@ -169,13 +169,17 @@ fhat = function(x_data, y_data, x, y_grid, p, q, mu, nu, h, kernel_type){
 
     # x constants
     x_scaled = sweep(x_sorted, 2, x)/(h^d)
-    if(check_inv(S_x(x_scaled, q, kernel_type)/(n*h^d))[1]== TRUE){
-      sx_mat = solve(S_x(x_scaled, q, kernel_type)/(n*h^d))
-    } else{
-      singular_flag = TRUE
-      sx_mat = matrix(0L, nrow = length(e_nu), ncol = length(e_nu))
+    if(length(x_scaled) == 0){
+      bx = 0
+    } else {
+      if(check_inv(S_x(x_scaled, q, kernel_type)/(n*h^d))[1]== TRUE){
+        sx_mat = solve(S_x(x_scaled, q, kernel_type)/(n*h^d))
+        } else{
+          singular_flag = TRUE
+          sx_mat = matrix(0L, nrow = length(e_nu), ncol = length(e_nu))
+          }
+      bx = b_x(x_scaled, sx_mat, e_nu, q, kernel_type)
     }
-    bx = b_x(x_scaled, sx_mat, e_nu, q, kernel_type)
 
     for (j in 1:ng){
       y = y_grid[j]
@@ -365,13 +369,18 @@ cov_hat = function(x_data, y_data, x, y_grid, p, q, mu, nu, h, kernel_type){
 
     # x constants
     x_scaled = sweep(x_sorted, 2, x)/(h^d)
-    if(check_inv(S_x(x_scaled, q, kernel_type)/(n*h^d))[1]== TRUE){
-      sx_mat = solve(S_x(x_scaled, q, kernel_type)/(n*h^d))
-    } else{
-      singular_flag = TRUE
-      sx_mat = matrix(0L, nrow = length(e_nu), ncol = length(e_nu))
+    if(length(x_scaled) == 0){
+      bx = 0
+    } else {
+      if(check_inv(S_x(x_scaled, q, kernel_type)/(n*h^d))[1]== TRUE){
+        sx_mat = solve(S_x(x_scaled, q, kernel_type)/(n*h^d))
+      } else{
+        singular_flag = TRUE
+        sx_mat = matrix(0L, nrow = length(e_nu), ncol = length(e_nu))
+      }
+      bx = b_x(x_scaled, sx_mat, e_nu, q, kernel_type)
     }
-    bx = b_x(as.matrix(x_scaled), sx_mat, e_nu, q, kernel_type)
+
 
     # initialize matrix
     c_hat = matrix(0L, nrow=ng, ncol=ng)

diff --git a/R/lpcde/man/lpbwcde.Rd b/R/lpcde/man/lpbwcde.Rd
diff --git a/R/lpcde_replication.R b/R/lpcde_replication.R
@@ -11,7 +11,7 @@ y_data = matrix(rnorm(n, mean=x_data, sd=1))
 y_grid = seq(from=-2, to=2, length.out=10)
 
 #Standard density estimation
-model1 = lpcde::lpcde(x_data=x_data, y_data=y_data, y_grid=y_grid, x=0, bw=1, rbc = TRUE)
+model1 = lpcde::lpcde(x_data=x_data, y_data=y_data, y_grid=y_grid, x=0, rbc = TRUE)
 summary(model1)
 #Bandwidth selection
 model2 = lpcde::lpbwcde(y_data=y_data, x_data=x_data, x=0, y_grid = y_grid)

diff --git a/todo.org b/todo.org
@@ -2,13 +2,13 @@
 #+author: Rajita Chandak
 
 * JOSS Revision
-** TODO bandwidth selection error example
+** DONE bandwidth selection error example
 *** DONE typo in documentation
 ** TODO check grid redundancy
 ** TODO adding predict method
-** TODO add more unit tests
-*** TODO checking integrates to 1
-*** TODO checking non-negativity
+** DONE add more unit tests
+*** DONE checking integrates to 1
+*** DONE checking non-negativity
 ** DONE ggplot deprecated warnings
 ** DONE multivariate warning
 ** DONE add missing doi