0

I want to sample theta[j] as first for j=1,2,...,71, then draw replicated( like 1000 times) yrep[k] form Bin(n[j], theta[j]), n[j] is known.

Then for each theta[j], I have 1000 yrep. If I want to choose the maximal element among this 1000 yrep. Finally, I will have 71 maximal elements.

My for loop in R is as follows:

theta<-NULL
yrep<-NULL
test<-NULL
k=1
for(i in 1:1000){
  for(j in 1:71){
    theta[j] <- rbeta(1,samp_A+y[j], samp_B+n[j]-y[j])
    yrep[k]<-rbinom(1, n[j], theta[j])
    k=k+1
  }
  t<-c(test, max(yrep))
}

I am not sure if this one is right.

 #Data
  y <- c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,
   2,1,5,2,5,3,2,7,7,3,3,2,9,10,4,4,4,4,4,4,4,10,4,4,4,5,11,12,
   5,5,6,5,6,6,6,6,16,15,15,9,4)
  n <- 
   c(20,20,20,20,20,20,20,19,19,19,19,18,18,17,20,20,20,20,19,19,18,18,25,24,
   23,20,20,20,20,20,20,10,49,19,46,27,17,49,47,20,20,13,48,50,20,20,20,20,
   20,20,20,48,19,19,19,22,46,49,20,20,23,19,22,20,20,20,52,46,47,24,14)


  #Evaluate densities in grid
  x <- seq(0.0001, 0.9999, length.out = 1000)


  #Compute the marginal posterior of alpha and beta in hierarchical model Use grid

  A <- seq(0.5, 15, length.out = 100)
  B <- seq(0.3, 45, length.out = 100)

  #Make vectors that contain all pairwise combinations of A and B

  cA <- rep(A, each = length(B))
  cB <- rep(B, length(A))

 #Use logarithms for numerical accuracy!

 lpfun <- function(a, b, y, n) log(a+b)*(-5/2) +
  sum(lgamma(a+b)-lgamma(a)-lgamma(b)+lgamma(a+y)+lgamma(b+n-y)- 
   lgamma(a+b+n))
lp <- mapply(lpfun, cA, cB, MoreArgs = list(y, n))

 #Subtract maximum value to avoid over/underflow in exponentiation

 df_marg <- data.frame(x = cA, y = cB, p = exp(lp - max(lp)))

 #Sample from the grid (with replacement)

  nsamp <- 100
  samp_indices <- sample(length(df_marg$p), size = nsamp,
                   replace = T, prob = df_marg$p/sum(df_marg$p))
  samp_A <- cA[samp_indices[1:nsamp]]
  samp_B <- cB[samp_indices[1:nsamp]]
   df_psamp <- mapply(function(a, b, x) dbeta(x, a, b),
               samp_A, samp_B, MoreArgs = list(x = x)) %>%
   as.data.frame() %>% cbind(x) %>% gather(ind, p, -x)
Hermi
  • 357
  • 1
  • 11

1 Answers1

0

You should provide us some data to work with, i didn't got to test it but i think this works:

for(j in 1:71){
  theta[j] <- rbeta(1,samp_A+y[j], samp_B+n[j]-y[j])
  for(i in 1:1000){
    yrep[i]<-rbinom(1, n[j], theta[j])}
  test = c(test, max(yrep))}
Ricardo Semião e Castro
  • 4,366
  • 1
  • 8
  • 27
  • @Lucas You're welcome! PS: you forgot adding data for `samp_A` and `samp_B`, and if my answer worked, please accept it to mark this post as solved. – Ricardo Semião e Castro Oct 16 '20 at 20:16
  • I test your code. I also test this code for the mean, min, and sample deviation. The p-value of max is 0.12 and the p-value of mean is 0.38. But why the p-values of min and sample deviation are 1 and 0? – Hermi Oct 16 '20 at 20:22
  • I also added the data for samp_A and samp_B. The p-value is just mean(test>=mean(y)) for mean, mean(test>=max(y)) for max, mean(test>=sd(y)) for sd. – Hermi Oct 16 '20 at 20:27
  • Without context of why you're calculating this statistics and what are those distributions i don't know how to interpret it, but it probably have something to do that P[y>=min(y)]=1 (a random variable is equal or higher than its minimal value with probability 1). If you have a question on the statistics, not on the code, it is best to ask this question in Cross Validated. PS: You said that the p-value for sd was 1 but for me it was `> mean(test>=sd(y)) [1] 0.9859155`. – Ricardo Semião e Castro Oct 16 '20 at 20:53
  • Thank you. Got it. – Hermi Oct 16 '20 at 23:26
  • Sorry, could you help me answer this similar question https://stackoverflow.com/questions/64399321/how-to-write-a-for-loop-to-compute-max-of-each-column-for-a-dataset-in-r? Thank you. – Hermi Oct 17 '20 at 05:34