[R] Repeating a R code and counting how many repetitions working

[varin sacha]

Dear Rui,

Many thanks, I got it now, thanks.

In my R code, I have moved the replication term 
N <-1000
out <- replicate(N, {

This time I get what I am looking for.

########################################
library(boot)

s<- sample(178:798, 10000, replace=TRUE)
mean(s)

N <- 10000
out <- replicate(N, {
a<- sample(s,size=5)
mean(a)
dat<-data.frame(a)
med<-function(d,i) {
temp<-d[i,,drop=TRUE]
mean(temp)
}

  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, 21:29:18 UTC+1, Rui Barradas <ruipbarradas using sapo.pt> a écrit : 





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 <ruipbarradas using 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 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.
>>