r - 如何在一个循环中选择最大元素在 r 中编写一个双循环?
问题描述
我想先对 j=1,2,...,71 采样 theta[j],然后从 Bin(n[j], theta[j]), n 中复制(比如 1000 次)yrep[k] [j] 是已知的。
然后对于每个 theta[j],我有 1000 yrep。如果我想在这 1000 yrep 中选择最大元素。最后,我将有 71 个最大元素。
我在 R 中的 for 循环如下:
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))
}
我不确定这个是否正确。
#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)
解决方案
你应该向我们提供一些数据来处理,我没有测试它,但我认为这可行:
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))}
推荐阅读
- python - 在两个 DataFrame 之间获得 delta 的最简单方法
- ruby - Ruby 方法有哪些形式?
- asp.net - DNN 模块弹出窗口显示在底部
- jquery - 影响所有用户的 Jquery 动态内容
- spring - STS 中缺少 Spring/Spring bean 配置文件选项
- shell - 使用 psql 命令行时引用字符串解释的问题
- ag-grid - 使用树数据进行行拖动
- javascript - 在 postgresql 数据库中存储 semver 版本字符串以进行范围查询
- laravel - Laravel - 您无权访问此服务器上的 /images/
- python - 无法正确识别实例列表中的 Python 类属性