高校数学の質問スレ(医者・東大卒専用) Part438 (991レス)
上下前次1-新
抽出解除 必死チェッカー(本家) (べ) 自ID レス栞 あぼーん
793: 132人目の素数さん [sage] 04/18(金)12:53 ID:XpMQ9r6R(1)
# Required package
library(HDInterval)
# Function to estimate alpha and beta of a Beta distribution
# from given HPI [L, U] and mean m
estimate_beta_params_from_HPI <- function(L, U, m, conf = 0.95, verbose = TRUE) {
# Initial guesses for alpha and beta
alpha0 <- m * 10
beta0 <- (1 - m) * 10
# Objective function: minimize error between estimated HPI/mean and given values
objective_fn <- function(par) {
a <- par[1]
b <- par[2]
if (a <= 0 || b <= 0) return(Inf)
hpi <- hdi(qbeta, shape1 = a, shape2 = b, credMass = conf)
est_mean <- a / (a + b)
hpi_error <- (hpi[1] - L)^2 + (hpi[2] - U)^2
mean_error <- (est_mean - m)^2
return(hpi_error + mean_error * 10) # Penalize deviation in mean
}
# Optimization
res <- optim(c(alpha0, beta0), objective_fn, method = "L-BFGS-B",
lower = c(0.001, 0.001))
alpha_hat <- res$par[1]
beta_hat <- res$par[2]
# Validate result
estimated_mean <- alpha_hat / (alpha_hat + beta_hat)
estimated_hpi <- hdi(qbeta, shape1 = alpha_hat, shape2 = beta_hat, credMass = conf)
if (verbose) {
cat("---- Result ----\n")
cat(sprintf("Estimated alpha: %.4f\n", alpha_hat))
cat(sprintf("Estimated beta : %.4f\n", beta_hat))
cat(sprintf("→ Mean : %.4f (target: %.4f)\n", estimated_mean, m))
cat(sprintf("→ %.0f%% HPI : [%.4f, %.4f] (target: [%.4f, %.4f])\n",
conf * 100, estimated_hpi[1], estimated_hpi[2], L, U))
}
return(list(alpha = alpha_hat,
beta = beta_hat,
mean = estimated_mean,
hpi = estimated_hpi))
}
# --- Example usage ---
# Suppose we are given:
# - Mean = 0.6
# - 95% HPI = [0.45, 0.75]
result <- estimate_beta_params_from_HPI(L = 0.45, U = 0.75, m = 0.6)
上下前次1-新書関写板覧索設栞歴
スレ情報 赤レス抽出 画像レス抽出 歴の未読スレ AAサムネイル
ぬこの手 ぬこTOP 1.568s*