Dear John,
Great, many thanks for your quick help.
Best,
SV
Le samedi 8 janvier 2022, 19:09:39 UTC+1, Fox, John <jfox at 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 at
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 at 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.