Вызов R2WinBUGS один раз для разных номеров, а не несколько раз для каждого номера
Я могу повторить байесовскую модель WinBUGS с n
лица, упорядоченные по s
, Для 50 человек, с последовательностью 10, я бы позвонил в WinBUGS 5 раз: для 10, 20, 30, 40 и 50 человек. Вместо того, чтобы запускать разные модели WinBUGS, я хочу использовать одну и ту же модель каждый раз. Я хочу использовать одну модель для 50 человек, а затем использовать эту модель для 10, 20, 30 и 40 человек. Другими словами, я хочу, чтобы R вызывал WinBUGS только один раз, а не пять раз.
Код, который я в настоящее время использую для 50 человек, с последовательностью 10:
n <- 50
s <- 10
bugs.model2.fn(n, s)
С выводом в виде графика, исследуем значения дельты:
Функция bugs.model2.fn(n,s)
описано ниже. Он выводит сюжет на основе применения sapply
с n
а также s
на другую функцию, которая генерирует данные с помощью WinBUGS (bugs.model1.fn(n)
)
library(R2WinBUGS)
library(ggplot2)
library(reshape2)
bugs.model1.fn <- function(n) {
m <- matrix(rep(0, 4), ncol=4) # dummy matrix
model <- "/Users/Model.txt"
# data
t1 <- rbinom(1, n, 0.95)
t2 <- rbinom(1, t1, 0.85)
t3 <- rbinom(1, t2, 0.95)
k1 <- n - t1
k2 <- t1 - t2
k3 <- t2 - t3
data <- list("k1", "k2", "k3", "t1", "t2", "t3")
# inits
myinits <- list(
list(theta1 = 0.1, theta2 = 0.1, theta3 = 0.1),
list(theta1 = 0.9, theta2 = 0.9, theta3 = 0.9))
parameters = c("theta1", "theta2", "theta3", "delta1", "delta2", "delta3")
# Bayesian model plot
bugs.model <- bugs(data, inits=myinits, parameters, model.file=model, bugs.directory=bugsdir,
n.chains=2, n.iter=10000, n.burnin=5000, n.thin=5, DIC=F, codaPkg=F, debug=F)
delta1 <- sum(bugs.model$sims.list$delta1 > 0.05) / length(bugs.model$sims.list$delta1)
delta2 <- sum(bugs.model$sims.list$delta2 > 0.05) / length(bugs.model$sims.list$delta2)
delta3 <- sum(bugs.model$sims.list$delta3 > 0.05) / length(bugs.model$sims.list$delta3)
r <- cbind(n, delta1, delta2, delta3)
mat <- rbind(m, r)
matrix <- mat[-1,]
return(matrix)
}
bugs.model2.fn <- function(n, s) {
output <- sapply(seq(from=max(2, s), to=n, by=s), FUN=function(n) return(bugs.model1.fn(n)))
df <- as.data.frame(t(output))
dfm <- melt(df, id.var = "n")
plot <- ggplot(dfm, aes(x = n, y = value, colour = variable)) + geom_point() + geom_line() +
labs(title = "% of delta difference >5%", x = "Number of individuals", y = "Proportion")
return(grid.arrange(plot))
}
Это использует Model.txt
, который содержит:
# Difference Between Two Rates
model {
# Prior on Rates
theta1 ~ dbeta(1,1)
theta2 ~ dbeta(1,1)
theta3 ~ dbeta(1,1)
# Observed Counts
k1 ~ dbin(theta1,t1)
k2 ~ dbin(theta2,t2)
k3 ~ dbin(theta3,t3)
# Difference between Rates
delta1 <- theta1-theta2
delta2 <- theta1-theta3
delta3 <- theta2-theta3
}
Я использую последнюю версию R2WinBUGS (версия 2.1 - 21), связанную с WinBUGS на MacOS High Sierra (версия 10.13.3) с Wine (версия 1.8.3) в качестве эмулятора Windows.