Bootstrap Confidence Intervals for Manhattan Distance

25 Views Asked by At

I would like to create 95% CI for the manhattan distance between each band for reflectance data that has been modified. The example bellow has 50 samples, which initial measurements will be stored as "or" and then it is modified version as "tr" in the city.fun. I wrote the following code but I am not sure if it is doing what I want it to do, which is bootstrap the value of each diagonal outcome (band similarity between before and after, per waveband (WL)).

library(spectrolab)
library(tidyverse)
library(boot)
library(flexclust)

spec = as_spectra(spec_matrix_example, name_idx = 1)
spec_mean = aggregate(spec, by = names(spec), mean, try_keep_txt(mean))
df = as.data.frame(spec, fix_names = 'none')
df_t = t(df)
colnames(df_t) <- df_t[1,]
df_t = df_t[-1,] 
df_2 = mutate_if(df, is.numeric, ~.*10)
df_2t = t(df_2)
colnames(df_2t) <- df_2t[1,]
df_2t = df_2t[-1,] 

sp_mx = data.frame(cbind(df_t, df_2t))
sp_mx = sp_mx %>% mutate_if(is.character, as.numeric)
sp_mx = as.matrix(sp_mx)

city.fun <- function(dat,idx) {
  or <- dat[idx,1:50]
  tr <- dat[idx,51:100]
  c(diag(dist2(or,tr, method = 'manhattan', p = 1)))
}

set.seed(2023)

bst <- boot(sp_mx, city.fun, R = 10)
bst
plot(bst)
boot.ci(boot.out = bst, type = c("norm", "basic","perc"))

I calculated the CI as the percentile CI

boot_stats <- data.frame(WL = as.numeric(rownames(data.frame(bst$t0))),
                         original = bst$t0, 
                         bias = colMeans(bst$t) - bst$t0, 
                         std_error = apply(bst$t,2,sd),
                         L.CI = apply(bst$t, 2, quantile,prob=0.025), 
                         U.CI = apply(bst$t, 2, quantile,prob=0.975)
)

ggplot(boot_stats) + 
  geom_point(aes(WL, original), size=0.5) + 
  geom_point(aes(WL,L.CI), size=0.5) + 
  geom_point(aes(WL,U.CI), size=0.5)

When I look at lower and upper limits it seems like the code is bootstrapping all the data together instead of for each band. Am I correct, or the code is doing what is intended to do? If the code is not doing so, how can I fix it?

0

There are 0 best solutions below