Hello,
As for the code, I think it's simple. Before calling boot.ci you need to
call boot. And replicate calls the expression between {} N times. After
running boot.ci, the component bca is extracted. See ?boot.ci, section
Value, after bca:
These latter four components will be matrices with 5 columns, the first
column containing the level, the next two containing the indices of the
order statistics used in the calculations and the final two the
calculated endpoints themselves.
The four components are all intervals but normal, which is a matrix with
3 columns only.
So I extract the 4th and 5th columns.
boot.ci(boot.out, type = "bca")$bca[, 4:5]
Also, try dropping the dimensions in the med function.
med <- function(d, i) {
temp <- d[i, , drop = TRUE]
mean(temp)
}
Your sample a is so small that there are only 3125 arrangements with
repetition, and only 111 different mean values of those arrangements.
arr <- RcppAlgos::permuteGeneral(a, length(a), repetition = TRUE)
dim(arr)
#[1] 3125 5
length(table( rowMeans(arr) ))
#[1] 111
Now, the equal values. The return value of 1 means that all BCa
intervals include mean(s). The outer mean(.) is a mean of an indicator
function given by a logical conjunction. And there are only TRUE's.
Hope this helps,
Rui Barradas
Às 15:37 de 11/12/21, varin sacha escreveu:
Dear Rui,
I really thank you a lot for your response. Even if I don't fully understand
your R code, there is something strange happening while running the code.
Indeed, I have run this R code here below 20 times and I always got the same
answer : [1] 1
########################################
library(boot)
s= sample(178:798, 10000, replace=TRUE)
mean(s)
a=sample(s,size=5)
mean(a)
dat<-data.frame(a)
med<-function(d,i) {
temp<-d[i,]
mean(temp)
}
N <- 100
out <- replicate(N, {
boot.out <- boot(data = dat, statistic = med, R = 10000)
boot.ci(boot.out, type = "bca")$bca[, 4:5]
})
mean(out[1,] < mean(s) & mean(s) < out[2,])
########################################
Le samedi 11 décembre 2021, 12:55:43 UTC+1, Rui Barradas <ruipbarra...@sapo.pt>
a écrit :
Hello,
The problem can be solved with ?replicate.
in the code below I only repeat N <- 1e3, not 1e4.
set.seed(2021)
N <- 1e3
out <- replicate(N, {
boot.out <- boot(data = dat, statistic = med, R = 10000)
boot.ci(boot.out, type = "bca")$bca[, 4:5]
})
mean(out[1,] < mean(s) & mean(s) < out[2,])
Hope this helps,
Rui Barradas
Às 10:47 de 11/12/21, varin sacha via R-help escreveu:
Dear R-experts,
Here below my R code. I am trying to do 2 things :
1) I would like to repeat this R code 10000 times
2) Out of the 10000 repetitions I would have liked R to tell me how many times the "true"
mean value of the population called "s" in my R example here below is included in the
10000 BCa confidence intervals constructed.
Many thanks for your precious help.
########################################
library(boot)
s= sample(178:798, 10000, replace=TRUE)
mean(s)
a=sample(s,size=5)
mean(a)
dat<-data.frame(a)
med<-function(d,i) {
temp<-d[i,]
mean(temp)
}
boot.out<-boot(data=dat,statistic=med,R=10000)
boot.ci(boot.out,type="bca")
########################################
______________________________________________
R-help@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.
______________________________________________
R-help@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.
