Вызов 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.

0 ответов

Другие вопросы по тегам