高校数学の質問スレ(医者・東大卒専用) Part438 (991レス)
上下前次1-新
抽出解除 必死チェッカー(本家) (べ) 自ID レス栞 あぼーん
リロード規制です。10分ほどで解除するので、他のブラウザへ避難してください。
815: 05/01(木)09:07 ID:L1qIlz9/(1/4) AAS
ディリクレ事前分布のパラメータαを階層化することで、より信頼性の高いベイズ推定が可能となる。
特にこの問題のように「実際に歪んでいる可能性がある」かつ「繰り返しが少ない」ケースでは、階層ベイズモデルはより適切な枠組みです。
816: 05/01(木)09:09 ID:L1qIlz9/(2/4) AAS
Stanで作ったらコンパイルに時間がかかる。簡単なモデルはJAGSの方がいい。離散変数も扱えるし。
# JAGS model
library(rjags)
# Prepare the data
outcome_data <- c(rep(1, 17), rep(2, 21), rep(3, 15), rep(4, 21), rep(5, 20), rep(6, 6))
N <- length(outcome_data)
data_jags <- list(outcome = outcome_data, N = N)
# Initial values (adjust as needed)
inits_jags <- list(
list(alpha = rep(1, 6), eta = 1),
list(alpha = runif(6, 0.1, 2), eta = 5)
)
# Compile the model
model_jags <- jags.model(
file = "hierarchical_dice_model.jag",
data = data_jags,
n.chains = 2,
n.adapt = 1000
)
# Sampling
samples_jags <- coda.samples(
model = model_jags,
variable.names = c("prob_simplex", "alpha", "eta"),
n.iter = 4000
)
# Summary of the results
cat("\nJAGS Sampling Results Summary:\n")
summary(samples_jags)
# Extract posterior samples (prob_simplex)
prob_simplex_posterior_jags <- as.matrix(samples_jags[, grep("prob_simplex", varnames(samples_jags))])
head(prob_simplex_posterior_jags)
# Plotting (example: posterior distribution of probabilities for each outcome)
cat("\nPosterior Distribution Plots for Each Outcome:\n")
par(mfrow = c(2, 3))
for (i in 1:6) {
plot(prob_simplex_posterior_jags[, i], type = "l", main = paste("Prob[", i, "]"), xlab = "Iteration", ylab = "Probability")
abline(h = 1/6, col = "red", lty = 2)
}
dice_prob_mean=prob_simplex_posterior_jags
colors <- c("skyblue", "lightcoral", "lightgreen", "gold", "lightsalmon", "lightcyan")
for (i in 1:ncol(dice_prob_mean)) {
BEST::plotPost(dice_prob_mean[, i],
compVal=1/6,
xlab=paste("pip ", i),
xlim=c(0, 0.4),
main="",
col=colors[i],
border="black")
}
817: 05/01(木)09:14 ID:L1qIlz9/(3/4) AAS
事後分布が出せたのであとはオッズ比などの計算も容易。
画像リンク
818: 05/01(木)11:18 ID:L1qIlz9/(4/4) AAS
ニッチな値の探索処理が終了しないコード
rm(list=ls())
library(fmsb)
library(parallel)
alpha <- 0.05
# Function to perform a single simulation
sim_single <- function(N = 100) {
A <- sample(1:N, 1)
while (A > N - 2) A <- sample(1:N, 1)
B <- sample(N - A - 1, 1)
C <- N - A - B
ABC <- c(A, B, C)
abc <- sapply(ABC, function(x) sample(x, 1))
x <- abc
n <- ABC
m <- rbind(s = x, f = n - x)
bonf_res <- pairwise.fisher.test(x, n, p.adj = 'bonf')
holm_res <- pairwise.fisher.test(x, n, p.adj = 'holm')
fdr_res <- pairwise.fisher.test(x, n, p.adj = 'fdr')
none_res <- pairwise.fisher.test(x, n, p.adj = 'none')
bonf <- min(bonf_res$p.value, na.rm = TRUE)
holm <- min(holm_res$p.value, na.rm = TRUE)
fdr <- min(fdr_res$p.value, na.rm = TRUE)
none <- min(none_res$p.value, na.rm = TRUE)
list(m = m, bonf = bonf, holm = holm, fdr = fdr, none = none)
}
# Function to find a result that meets the criteria using parallel processing
find_result_parallel_loop <- function(alpha = 0.05, num_cores = detectCores() - 1) {
# Create a cluster of worker processes
cl <- makeCluster(num_cores)
# Ensure the cluster is stopped when the function exits
on.exit(stopCluster(cl))
# Export the sim_single function to the worker processes
clusterExport(cl, "sim_single")
# Load the fmsb library in the worker processes
clusterEvalQ(cl, library(fmsb))
cat("Searching for a result that meets the criteria...\n")
while (TRUE) {
# Run simulations in parallel
results <- parLapply(cl, 1:num_cores, function(i) { # Run as many simulations as cores
sim_single()
})
# Check the results for the desired condition
for (res in results) {
if (res$bonf > alpha && res$holm < alpha) {
cat("Result meeting the criteria found:\n")
return(res)
}
}
# If no result found in this batch, continue the loop
}
}
# Find the result using parallel processing until found
res_parallel_loop <- find_result_parallel_loop(alpha = alpha)
# Output the result (will be printed within the loop when found)
print(res_parallel_loop)
上下前次1-新書関写板覧索設栞歴
スレ情報 赤レス抽出 画像レス抽出 歴の未読スレ
ぬこの手 ぬこTOP 0.043s