[R] Percentile bootstrap for the median : error message

[varin sacha]

Dear John,

Great, many thanks for your quick help.

Best,
SV



Le samedi 8 janvier 2022, 19:09:39 UTC+1, Fox, John <jfox using mcmaster.ca> a écrit : 





Dear Sacha,

Here's your corrected and cleaned-up code:

> library(boot)
> set.seed(123)
> s <- rnorm(100000,0,1)
> (m <- median(s)) 
[1] 0.000946463

> med <- function(d,i) {
+  median(d[i, ])
+ }

> set.seed(456)
> N <- 100
> n<-5

> out <- replicate(N, {
+  dat <- data.frame(sample(s,size=n))
+  boot.out <- boot(data = dat, statistic = med, R = 10000)
+  boot.ci(boot.out, type = "perc")$perc[, 4:5]  
+ })

> mean(out[1, ] < m & m < out[2, ])
[1] 0.94

A couple of comments:

(1) I moved the definition of med() outside of the call to replicate() so the function doesn't get defined repeatedly -- something that, if I'm not mistaken, you were advised to do the last time you asked a very similar question.

(2) Maybe it's time to polish your debugging skills. I quickly found the error in the call to boot.ci() by calling browser() before the line dat <- data.frame(sample(s,size=n)) and stepping through the commands.

I hope this helps,
John

  -----------------------------
  John Fox, Professor Emeritus
  McMaster University
  Hamilton, Ontario, Canada
  Web: http::/socserv.mcmaster.ca/jfox

> On Jan 8, 2022, at 12:04 PM, varin sacha via R-help <r-help using r-project.org> wrote:
> 
> Dear R-experts,
> 
> Here below my R code for the percentile bootstrap confidence intervals with an error message. 
> Is there a way to make my R code work ?
> Many thanks for your help and time.
> 
> ############################################
> library(boot)
> 
> s=rnorm(100000,0,1)
> (m<-median(s)) 
> 
> N <- 100
> n<-5
> out <- replicate(N, {
> 
> dat<-data.frame(sample(s,size=n))
> med<-function(d,i) {
> median(d[i, ])
> }
> 
>  boot.out <- boot(data = dat, statistic = med, R = 10000)
> 
>  boot.ci(boot.out, type = "per")$per[, 4:5]
> 
> })
> 
> mean(out[1, ] < m & m < out[2, ])
> ############################################
> 
> ______________________________________________
> R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.