[過去ﾛｸﾞ] 高校数学の質問スレ（医者・東大卒専用） Part438 (1002ﾚｽ)
上下前次1-新
このｽﾚｯﾄﾞは過去ﾛｸﾞ倉庫に格納されています｡
次ｽﾚ検索 歴削→次ｽﾚ 栞削→次ｽﾚ 過去ﾛｸﾞﾒﾆｭｰ
839: 05/17(土)07:49 ID:Vpav5/5q(1) AAS
library(meta)
library(gemtc)
# データフレームの作成
data <- data.frame(
study = factor(c(1, 1, 2, 2)),
treatment = factor(c("Placebo", "Old", "Old", "New")),
events = c(19, 5, 5, 9),
total = c(202, 201, 204, 203)
)
# 治療効果の比較リスト
comparisons <- c("New - Placebo")
# ネットワークオブジェクトの作成
network <- mtc.network(data = data)
# ランダム効果モデルによるメタアナリシス
model <- mtc.model(network, type = "consistency", likelihood = "binom", link = "logit", linearModel = "random")
# モデルの実行
results <- mtc.run(model)
# 結果の概要
summary(results)
# 新薬 vs プラセボ の効果量（オッズ比）と95%信頼区間
relative.effect(results, "New", "Placebo")
# 新薬がプラセボより劣る確率（P(OR > 1)）を推定
prob.superiority <- 1 - pnorm(0, mean = results$sol[grep("d.New", names(results$sol))], sd = results$se[grep("d.New", names(results$se))])
cat("新薬がプラセボより劣る確率 (間接推定):", prob.superiority, "\n")
# 新薬がプラセボに対して有意に有効かどうかの評価
# オッズ比の95%信頼区間に1が含まれるかどうかで判断
or.new.vs.placebo <- exp(results$sol[grep("d.New", names(results$sol))])
ci.lower <- exp(results$sol[grep("d.New", names(results$sol))] - 1.96 * results$se[grep("d.New", names(results$se))])
ci.upper <- exp(results$sol[grep("d.New", names(results$sol))] + 1.96 * results$se[grep("d.New", names(results$se))])
cat("新薬 vs プラセボ オッズ比 (間接推定):", or.new.vs.placebo, "\n")
cat("95%信頼区間:", ci.lower, "-", ci.upper, "\n")
if (ci.upper < 1) {
cat("新薬はプラセボに対して有意に有効である可能性が高いです。\n")
} else if (ci.lower > 1) {
cat("新薬はプラセボに対して有意に劣る可能性があります。\n")
} else {
cat("新薬とプラセボの間に有意な差は見られない可能性があります。\n")
}
840: 05/20(火)23:31 ID:gwaBTE4C(1) AAS
library(R2jags)
# データ
data <- list(
nA1 = 100, rA1 = 80, # Study1: 治療A
nB1 = 100, rB1 = 40, # Study1: 治療B
nA2 = 100, rA2 = 10, # Study2: 治療A
nC2 = 100, rC2 = 5 # Study2: 治療C
)
# JAGSモデル（textConnection使用）
model_code <- "
model {
# 尤度関数
rA1 ~ dbin(pA1, nA1)
rB1 ~ dbin(pB1, nB1)
rA2 ~ dbin(pA2, nA2)
rC2 ~ dbin(pC2, nC2)
# 治療Aの階層モデル
mu_A ~ dbeta(1, 1)
tau_A ~ dgamma(0.001, 0.001)
pA1 ~ dbeta(mu_A * tau_A, (1 - mu_A) * tau_A)
pA2 ~ dbeta(mu_A * tau_A, (1 - mu_A) * tau_A)
sigma_A <- 1 / sqrt(tau_A) # SDに変換
# 治療BとCも階層化（平均リスクを別々に推定）
mu_B ~ dbeta(1, 1)
mu_C ~ dbeta(1, 1)
pB1 ~ dbeta(mu_B * 100, (1 - mu_B) * 100) # 高い精度を仮定
pC2 ~ dbeta(mu_C * 100, (1 - mu_C) * 100)
# リスク差
RD_A1_B1 <- pA1 - pB1
RD_A2_C2 <- pA2 - pC2
RD_B1_C2 <- pB1 - pC2 # B vs Cの直接比較
}
"
# JAGS実行
jags_model <-
(textConnection(model_code),
data = data, n.chains = 3, quiet=TRUE)
update(jags_model, 3000) #, progress.bar="none")
jags_samples <- coda.samples(jags_model,
c("mu_A", "sigma_A", "RD_A1_B1", "RD_A2_C2", "RD_B1_C2", "pA1", "pA2", "pB1", "pC2"),
n.iter=10000) # , progress.bar="none")
gelman.plot(jags_samples)
plot(jags_samples)
summary(jags_samples)
jags_samples |> as.matrix() |> as.data.frame() -> js
names(js)
841: 05/23(金)10:58 ID:or+7Cxzr(1) AAS
#
"
Construct a Monte Carlo study that investigates how the probability of coverage depends on the sample size and true proportion value. In the study, let n be 10, 25, 50, and 100 and let p be .05, .25, and .50. Write an R function that has three inputs, n, p, and the number of Monte Carlo simulations m,and will output the estimate of the exact coverage probability.
Implement your function using each combination of n and p and m = 1000 simulations.
Describe how the actual probability of coverage of the traditional interval depends on the sample size and true proportion value.
"
f = \(n,p,m=1000){
y=rbinom(m,n,p)
phat=y/n
se=sqrt(phat*(1-phat)/n)
lo=phat - qnorm(0.975)*se
up=phat + qnorm(0.975)*se
mean(lo < p & p < up)
}
f=Vectorize(f)
n_values = c(10, 25, 50,100)
p_values = c(0.05, 0.25, 0.5)
set.seed(123)
outer(n_values,p_values,f)
842: 05/24(土)02:35 ID:VetM3rz7(1/5) AAS
LearnBayes::beta.selectをoptimを使って算出
beta.optim <- function(x1, p1, x2, p2, verbose = TRUE) {
# -------------------------
# モーメント近似による初期値推定
# -------------------------
mu0 <- (x1 + x2) / 2 # 仮の平均
sigma2 <- ((x2 - x1) / 4)^2 # 仮の分散（中間50%幅から）
v <- mu0 * (1 - mu0) / sigma2 - 1
a0 <- mu0 * v
b0 <- (1 - mu0) * v
start <- c(a0, b0)
# -------------------------
# 最適化対象の誤差関数定義
# -------------------------
objective <- function(params) {
a <- params[1]
b <- params[2]
# pbeta による累積確率との差を2乗誤差で評価
err1 <- pbeta(x1, a, b) - p1
err2 <- pbeta(x2, a, b) - p2
return(err1^2 + err2^2)
}
# -------------------------
# 最適化（境界付き）
# -------------------------
result <- optim(
par = start,
fn = objective,
method = "L-BFGS-B",
lower = c(1e-4, 1e-4) # ベータ分布は a, b > 0
)
# -------------------------
# 結果の取り出し
# -------------------------
a_hat <- result$par[1]
b_hat <- result$par[2]
if (verbose) {
cat("推定されたパラメータ: a =", round(a_hat, 4),
", b =", round(b_hat, 4), "\n")
cat("pbeta(x1) =", round(pbeta(x1, a_hat, b_hat), 4),
"（目標:", p1, "）\n")
cat("pbeta(x2) =", round(pbeta(x2, a_hat, b_hat), 4),
"（目標:", p2, "）\n")
}
# -------------------------
# 結果を返す
# -------------------------
return(list(a = a_hat, b = b_hat, ss_value = result$value))
}
843: 05/24(土)08:17 ID:VetM3rz7(2/5) AAS
library(rjags)
# Fit a Bayesian logistic regression model using JAGS and return predictions and posterior summaries
fit_bayesian_logistic_jags <- function(data, formula, newdata,
n.chains = 3, n.iter = 5000, n.burnin = 1000) {
# Extract response variable name from the formula
response_var <- all.vars(formula)[1]
y <- data[[response_var]]
# Convert factor response to binary numeric (0/1)
if (is.factor(y)) y <- as.numeric(y) - 1
y <- as.numeric(y)
# Construct design matrices for training and new data
X <- model.matrix(formula, data)
new_X <- model.matrix(delete.response(terms(formula)), newdata)
# Prepare data list for JAGS
jags_data <- list(
y = y,
X = X,
n = nrow(X),
p = ncol(X),
new_X = new_X,
scale_beta = rep(2.5, ncol(X)) # Prior scale for each coefficient
)
# Define the JAGS model
model_string <- "
model {
for (j in 1:p) {
beta[j] ~ dt(0, 1 / pow(scale_beta[j], 2), 1)
}
for (i in 1:n) {
logit_p[i] <- inprod(X[i,], beta[])
y[i] ~ dbern(1 / (1 + exp(-logit_p[i])))
}
new_logit <- inprod(new_X[1,], beta[])
new_p <- 1 / (1 + exp(-new_logit))
}
"
# Initialize and run the JAGS model
model <- jags.model(textConnection(model_string), data = jags_data,
n.chains = n.chains, quiet = TRUE)
update(model, n.burnin)
# Draw posterior samples
samples <- coda.samples(model, c("beta", "new_p"), n.iter - n.burnin)
mat <- as.matrix(samples)
# Return results
list(
model = samples,
predicted_prob = mean(mat[, "new_p"]),
summary = summary(samples)
)
}
844: 05/24(土)08:18 ID:VetM3rz7(3/5) AAS
# Example data
data <- data.frame(
donation = c(0, 1000, 2000, 0, 3000, 0, 4000, 0, 5000, 0),
score = c(90, 40, 35, 88, 30, 85, 25, 92, 20, 89),
parent = c(0, 1, 1, 0, 1, 0, 1, 0, 1, 0),
admission = as.factor(c(0, 1, 1, 0, 1, 0, 1, 0, 1, 0))
)
# New observation to predict
newdata <- data.frame(
donation = 2500,
score = 40,
parent = 1
)
# Fit model and obtain results
set.seed(123)
result <- fit_bayesian_logistic_jags(
data = data,
formula = admission ~ donation + score + parent,
newdata = newdata
)
# Extract variable names including intercept
var_names <- colnames(model.matrix(admission ~ donation + score + parent, data))
# Extract beta coefficient summaries
beta_stats <- result$summary$statistics[grep("^beta\\[", rownames(result$summary$statistics)), c("Mean", "SD")]
beta_quants <- result$summary$quantiles[grep("^beta\\[", rownames(result$summary$quantiles)), c("2.5%", "97.5%")]
# Rename row names using variable names
rownames(beta_stats) <- var_names
rownames(beta_quants) <- var_names
# Display results
print(beta_stats)
print(beta_quants)
cat("Predicted probability:", round(result$predicted_prob, 3), "\n")
845: 05/24(土)08:37 ID:VetM3rz7(4/5) AAS
library(rjags)
# Fit a Bayesian logistic regression model using JAGS and return predictions and posterior summaries
fit_bayesian_logistic_jags <- function(data, formula, newdata,
n.chains = 3, n.iter = 5000, n.burnin = 1000) {
# Extract response variable name from the formula
response_var <- all.vars(formula)[1]
y <- data[[response_var]]
# Convert factor response to binary numeric (0/1)
if (is.factor(y)) y <- as.numeric(y) - 1
y <- as.numeric(y)
# Construct design matrices for training and new data
X <- model.matrix(formula, data)
new_X <- model.matrix(delete.response(terms(formula)), newdata)
# Prepare data list for JAGS
jags_data <- list(
y = y,
X = X,
n = nrow(X),
p = ncol(X),
new_X = new_X,
scale_beta = rep(2.5, ncol(X)) # Prior scale for each coefficient
)
# Define the JAGS model
model_string <- "
model {
for (j in 1:p) {
beta[j] ~ dt(0, 1 / pow(scale_beta[j], 2), 1)
}
for (i in 1:n) {
logit_p[i] <- inprod(X[i,], beta[])
y[i] ~ dbern(1 / (1 + exp(-logit_p[i])))
}
new_logit <- inprod(new_X[1,], beta[])
new_p <- 1 / (1 + exp(-new_logit))
}
"
# Initialize and run the JAGS model
model <- jags.model(textConnection(model_string), data = jags_data,
n.chains = n.chains, quiet = TRUE)
update(model, n.burnin)
# Draw posterior samples
samples <- coda.samples(model, c("beta", "new_p"), n.iter - n.burnin)
mat <- as.matrix(samples)
# Return results
list(
model = samples,
predicted_prob = mean(mat[, "new_p"]),
summary = summary(samples)
)
}
846: 05/24(土)21:16 ID:VetM3rz7(5/5) AAS
# dbeta(L,a,b) == dbbeta(U,a,b)
# Solve[L^(a-1)(1-L)^(b-1)==U^(a-1)(1-U)^(b-1), b]
L=1/7
U=1/5
credMass = 0.95
f = function(a) 1 + ((a - 1) * log(U / L)) / log((1 - L) / (1 - U))
g = function(a) pbeta(U,a,f(a)) - pbeta(L,a,f(a)) - credMass
(re=uniroot(g,c(1,1e5)))
curve(g(x),1,150,bty="l") ; abline(h=0,lty=3)
c(re$root,f(re$root))
847(1): 05/25(日)04:22 ID:P4nhnL8B(1/6) AAS
# dbeta(L,a,b) == dbbeta(U,a,b)
# Solve[L^(a-1)(1-L)^(b-1)==U^(a-1)(1-U)^(b-1), b]
L=1/7
U=1/5
credMass = 0.95
f = function(a) 1 + ((a - 1) * log(U / L)) / log((1 - L) / (1 - U))
g = function(a) pbeta(U,a,f(a)) - pbeta(L,a,f(a)) - credMass
(re=uniroot(g,c(1,1e5)))
a=re$root
b=f(a)
c(a,b)
curve(g(x),1,1.5*a,bty="l") ; abline(h=0,lty=3)
a/(a+b) # mean
(a-1)/(a-1+b-1) # mode
library(LearnBayes)
beta.select(list(x=1/7,p=0.025),list(x=1/5,p=0.975))
848: 05/25(日)05:57 ID:P4nhnL8B(2/6) AAS
＞847のRのコードをChatGPTで
Mathematicaにコメント付きで移植
(*
betaParameter 関数：
指定された信頼区間 [L, U] に、指定された信頼度 credMass（例: 95%）の確率質量を持つ
ベータ分布のパラメータ α, β を算出する。
*)
betaParameter[L_: 1/7, U_: 1/5, credMass_: 0.95] := Module[
{α, β},
(*
f[α] は、PDF[BetaDistribution[α, β], L] == PDF[BetaDistribution[α, β], U]
を満たすように β を α に基づいて算出する関数。
*)
f[α_] := 1 + ((α - 1) * Log[U / L]) / Log[(1 - L) / (1 - U)];
(*
g[α] は、ベータ分布 Beta[α, f[α]] の区間 [L, U] に
含まれる確率（CDFの差）を返す関数。
*)
g[α_] := CDF[BetaDistribution[α, f[α]], U] - CDF[BetaDistribution[α, f[α]], L];
(* g[α] = credMass を満たす α を数値的に求める *)
α = α /. FindRoot[g[α] == credMass, {α, 1, 1*^5}];
(* 対応する β を算出 *)
β = f[α];
(* 結果を返す *)
{α, β}
]
(* 関数を実行して α, β を取得 *)
{α, β} = betaParameter[]
(* g[α] を評価して、[L, U] に credMass の質量があることを確認 *)
g[α]
849: 05/25(日)06:42 ID:P4nhnL8B(3/6) AAS
>>847
このプロトタイプをAIに与えて描画機能やコメントをつけてもらった。
beta.parameter <- function(lower, upper, credMass = 0.95, verbose = FALSE) {
# Helper function to convert decimal numbers to fraction strings using MASS::fractions
fractionStr <- function(x) {
as.character(MASS::fractions(x))
}
# Function to compute beta parameter (beta) based on alpha,
# derived from the condition on the shape of the distribution between lower and upper
f <- function(alpha) {
1 + ((alpha - 1) * log(upper / lower)) / log((1 - lower) / (1 - upper))
}
# Root-finding function: difference between desired credible mass
# and the Beta CDF probability between lower and upper with parameters (alpha, f(alpha))
g <- function(alpha) {
pbeta(upper, alpha, f(alpha)) - pbeta(lower, alpha, f(alpha)) - credMass
}
# Find the root of g(alpha) = 0 over the interval [1, 1e5]
# to find the alpha value that satisfies the credible mass condition
re <- uniroot(g, c(1, 1e5))
alpha <- re$root
beta <- f(alpha)
# Calculate the mean of the Beta distribution
mean <- alpha / (alpha + beta)
# Calculate the mode if defined (alpha > 1 and beta > 1),
# otherwise set mode to NA as it is undefined
mode <- if (alpha > 1 && beta > 1) {
(alpha - 1) / (alpha + beta - 2)
} else {
NA
}
850: 05/25(日)06:42 ID:P4nhnL8B(4/6) AAS
# If verbose flag is TRUE, plot the Beta distribution and annotate results
if (verbose) {
# Generate x values from 0 to 1 for plotting the density
x <- seq(0, 1, length.out = 1000)
# Compute Beta density values at x
y <- dbeta(x, alpha, beta)
# Color bars within the credible interval [lower, upper] as "lightcoral",
# others as light gray ("gray70")
col <- ifelse(x >= lower & x <= upper, "lightcoral", "gray70")
# Plot histogram-like vertical lines representing the Beta density
plot(x, y, type = "h", col = col, lwd = 2,
main = sprintf("Beta(%.2f, %.2f) Distribution\n[%s, %s] with %.0f%% Probability Mass",
alpha, beta, fractionStr(lower), fractionStr(upper), credMass * 100),
xlab = "x", ylab = "Density", bty = "n")
# Add vertical dashed line for the mean, colored skyblue
abline(v = mean, col = "skyblue", lwd = 1, lty = 2)
# Add vertical dashed line for the mode if defined, colored dark green
if (!is.na(mode)) {
abline(v = mode, col = "darkgreen", lwd = 1, lty = 2)
}
# Prepare legend labels for mean, mode (if exists), and credible interval
labels <- c(
paste0("Mean = ", round(mean, 3)),
if (!is.na(mode)) paste0("Mode = ", round(mode, 3)) else NULL,
paste0("95% Credible Interval [", fractionStr(lower), ", ", fractionStr(upper), "]")
)
# Corresponding colors for legend items
colors <- c(
"skyblue",
if (!is.na(mode)) "darkgreen" else NULL,
"lightcoral"
)
# Line types for legend items (dashed lines for mean and mode, none for credible interval)
ltys <- c(2, if (!is.na(mode)) 2 else NULL, NA)
# Plot characters for legend (none for lines, solid square for interval)
pchs <- c(NA, if (!is.na(mode)) NA else NULL, 15)
# Add legend at the top of the plot with no box, scaled text size
legend("top", legend = labels, col = colors, bty = "n", cex = 0.9,
lty = ltys, pch = pchs)
}
# Return a named vector of calculated parameters and statistics
c(alpha = alpha, beta = beta, mean = mean, mode = mode)
}
851: 05/25(日)06:51 ID:P4nhnL8B(5/6) AAS
今回スムーズに機能拡張できた理由は、以下のようにプロトタイプが非常に明快かつ健全だったことが大きな要因です。
✅ プロトタイプの良さが拡張性を支えた理由
1. 構造がシンプル
中心となる数式（pbeta(U, a, b) - pbeta(L, a, b) = credMass）が明確で、ロジックが一貫していました。
f(a) で b を a の関数として定義しており、探索空間を1次元に抑えていた点も効率的。
2. 関数分離と再利用が可能
f(a) や g(a) が関数として定義されていたので、視覚化やバリエーションの追加が簡単でした。
3. 出力が明確
c(a, b) と返す構造が明瞭で、追加情報（期待値・最頻値など）も素直に足せました。
4. 数値的安定性が良好
uniroot() による根の探索は数値計算としても信頼性が高く、実用的。
852: 05/25(日)10:38 ID:P4nhnL8B(6/6) AAS
HDI_discrete <- function(prob_data, credMass) {
x = prob_data[, 1]
p = prob_data[, 2]
n = length(x)
sp = sort(p, index.return = TRUE)
i = sp$ix[seq(n, 1, -1)]
ps = p[i]
xs = x[i]
cp = cumsum(ps)
ii = 1:n
j = ii[cp >= credMass]
j = j[1]
achieved_credMass = cp[j]
hdi_set = sort(xs[1:j])
v = list(credMass = achieved_credMass, set = hdi_set)
return(v)
}
853: 05/29(木)17:36 ID:WPcwJ6cn(1) AAS
"
サイコロを50回なげて4回１の目がでた。
１の目のでる確率は1/6である仮説が正しい確率を求めよ。
計算に必要なモデルは適宜設定してよい。
例：サイコロを投げる前のこの仮説が正しい確率は一様分布に従う。
a/(a+b)=p
a=p(a+b)
(1-p)a=pb
b=a*(1-p)/p
a=b*p/(1-p)
a/b=p/(1-p)
"
rm(list=ls())
p=1/6
n=50
y=4
pbinom(y,n,p)
dbinom(y,n,p)
fn <- function(prob){
sub <-function(alpha) pbetat(p,prob,c(1,(1-p)/p)*alpha,c(y,n-y))$post
optimize(sub,c(1,100))$objective
}
fn = Vectorize(fn)
fn(0.5)
prior.probs=runif(1e6)
post.probs=fn(prior.probs)
summary(post.probs)
HDInterval::hdi(post.probs)
source("plotpost.R")
plotpost(post.probs,col="lightcoral")
lines(density(post.probs))
854: 05/30(金)06:24 ID:zoAuXcvc(1/3) AAS
par(bty="l")
y = c(43, 24, 100, 35, 85)
yn = max(y)
n = length(y)
B = 200
like = numeric(B)
for(i in yn:B) { like[i] = 1 / (i^n) }
pmf = like / sum(like)
mean = sum((1:B) * pmf)
plot(pmf, main="Posterior PMF of N", xlab="N", ylab="Probability", type="h", col="blue")
plot(cumsum(pmf), main="Posterior CDF of N", xlab="N", ylab="Cumulative Probability", type="s", col="blue")
abline(h=0.95, lty=3)
c(lower=yn, mean=mean, upper=which(cumsum(pmf) > 0.95)[1])
sd = sqrt(sum(((1:B) - mean)^2 * pmf))
prob_N_gt_150 = sum(pmf[151:B])
cat("Posterior Mean:", mean, "\n")
cat("Posterior Standard Deviation:", sd, "\n")
cat("P(N > 150):", prob_N_gt_150, "\n")
855: 05/30(金)19:44 ID:zoAuXcvc(2/3) AAS
(* pbetat関数の定義 *)
pbetat[p0_, prob_, ab_, data_] := Module[{a, b, s, f, lbf, bf, post},
a = ab[[1]];
b = ab[[2]];
s = data[[1]];
f = data[[2]];
lbf = s * Log[p0] + f * Log[1 - p0] + Log@Beta[a, b] - Log@Beta[a + s, b + f];
bf = Exp[lbf];
post = prob * bf / (prob * bf + 1 - prob);
<|"bf" -> bf, "post" -> post|>
]
(* 関数fの定義 *)
f[alpha_] := pbetat[1/6, 0.5, {alpha, 5*alpha - 4}, {4, 50 - 4}]["post"]
(* 最小化 *)
result = NMinimize[{f[alpha], 0 <= alpha <= 20}, alpha]
856: 05/30(金)20:13 ID:zoAuXcvc(3/3) AAS
p_post_null <- function(p0, prior, alpha, beta, success, failure){
# Calculate the total number of trials from successes and failures.
total = success + failure
# Calculate the likelihood of the data under the null hypothesis (H0).
# This assumes a binomial distribution where the success probability is p0.
m0 = dbinom(success, total, p0)
# Calculate the marginal likelihood of the data under the alternative hypothesis (H1).
# Under H1, the success probability is assumed to follow a Beta distribution
# with parameters alpha and beta. This function (from the 'extraDistr' package)
# computes the marginal likelihood by integrating over all possible success probabilities.
m1 = extraDistr::dbbinom(success, total, alpha, beta)
# Calculate the Bayes Factor (BF01).
# This is the ratio of the likelihood under the null hypothesis (m0)
# to the marginal likelihood under the alternative hypothesis (m1).
BF01 = m0 / m1
# Calculate the posterior probability of the null hypothesis.
# This updates the prior belief (prior) based on the evidence from the data (BF01).
p_post = prior * BF01 / (prior * BF01 + 1 - prior)
# Return both the Bayes Factor and the posterior probability as a named vector.
c(BF01 = BF01, p_post = p_post)
}
# Optimize the 'alpha' parameter to maximize the posterior probability of the null hypothesis.
# We're trying to find the 'alpha' value (within the range 0 to 1e6) that makes the
# null hypothesis most plausible, given the data and the relationship beta = 5*alpha - 4.
# p_post_null(...)[2] specifically extracts the 'p_post' value from the function's output.
optimize(function(alpha) p_post_null(p0 = 1/6,
prior = 1/2,
alpha,
5 * alpha - 4, # Beta parameter is a function of alpha
success = 4,
failure = 50 - 4)[2],
c(0, 1e6)) # Search range for alpha
# Calculate the Bayes Factor and posterior probability using specific alpha and beta values.
# These specific values (alpha = 50/9, beta = 5*50/9 - 4) are likely the result of the
# optimization step above, or pre-determined values that are of interest.
p_post_null(1/6, 1/2, 50/9, 5 * 50/9 - 4, 4, 50 - 4)
857: 05/31(土)05:25 ID:jzcOJBMt(1) AAS
#' @title ベイズ事後確率計算関数
#' @description 帰無仮説と対立仮説の事後確率を計算
#' @param s 観測成功数（1の目が出た回数）
#' @param n 総試行回数
#' @param p0 帰無仮説の確率（例: 1/6）
#' @param prior 帰無仮説の事前確率（0～1）
#' @param alpha 対立仮説のベータ分布αパラメータ
#' @param beta 対立仮説のベータ分布βパラメータ
#' @return list(bf01=ベイズファクター, post=事後確率, method=使用手法)
calculate_posterior <- function(s, n, p0, prior, alpha, beta) {
# 入力検証
stopifnot(
s >= 0, n > 0,
p0 > 0 && p0 < 1,
prior >= 0 && prior <= 1,
alpha > 0, beta > 0
)
# 帰無仮説の尤度計算
m0 <- dbinom(s, n, p0)
# 対立仮説の周辺尤度計算（extraDistr有無で自動切替）
if (requireNamespace("extraDistr", quietly = TRUE)) {
m1 <- extraDistr::dbbinom(s, n, alpha, beta)
method <- "extraDistr::dbbinom()"
} else {
integrand <- function(p) dbinom(s, n, p) * dbeta(p, alpha, beta)
m1 <- integrate(integrand, 0, 1)$value
method <- "数値積分"
}
# ベイズファクターと事後確率計算（指定された式を使用）
bf01 <- m0 / m1
post <- m0 * prior / (m0 * prior + (1 - prior) * m1)
# 結果をリストで返す
list(
bf01 = bf01,
post = post,
method = method,
inputs = list(s = s, n = n, p0 = p0, prior = prior, alpha = alpha, beta = beta),
likelihoods = list(m0 = m0, m1 = m1)
)
}
858: 06/01(日)07:35 ID:/PZFDI/g(1/2) AAS
auc = integrate(unnorm_posterior,0,Inf)$value
integrate(\(x) x*unnorm_posterior(x)/auc, 0,Inf)
log_auc <- log(integrate(\(x) exp(log(unnorm_posterior(x))), 0, Inf)$value)
integrate(\(x) x * exp(log(unnorm_posterior(x)) - log_auc), 0, Inf)$value
859: 06/01(日)09:23 ID:9o1m2vAK(1) AAS
Calculate the marginal likelihood of the data under the alternative hypothesis (H1).
wwwwwwwwwwwwwwwwwwwwwwwwwwwwww
あほ～～
wwwwwwwwwwwwwwwwwwwwwwwwwwwwww
860: 06/01(日)10:41 ID:/PZFDI/g(2/2) AAS
はい、Gelmanらの推奨は確かに進化しており、近年ではハーフt分布（half-t）がより一般的に推奨されています。この変化の背景と具体的な推奨内容を、理論的根拠と実践的なアドバイスに分けて解説します。
1. Gelmanの推奨の変遷
(1) 初期（2006年頃）: ハーフコーシー（自由度ν=1）の推奨
推奨理由:
分散パラメータの弱情報事前分布として適切
裾が重く、極端な値への過剰なペナルティを回避
問題点:
自由度ν=1ではスケールの選択に敏感
小規模データで不安定になる可能性
(2) 近年（2013年以降）: 自由度ν>1のハーフt分布へ移行
推奨理由:
自由度ν=3〜7でより安定した推定が可能
スケールパラメータの影響が緩和
代表的な論文:
Gelman et al. (2013), "The Prior Can Often Only Be Understood in the Context of the Likelihood"
861: 06/02(月)11:14 ID:GMuHFUYr(1/2) AAS
x = c(-0.86, -0.3, -0.05, 0.73)
n = c(5, 5, 5, 5)
y = c(0, 1, 3, 5)
(data = cbind(x, n, y))
(response = cbind(y, n - y) )
results = glm(response ~ x, family = binomial)
#summary(results)
-results$coef[1]/results$coef[2]
library(MASS) # mvrnorm を使うため
# 推定された係数と共分散行列
beta_hat = coef(results)
(vcov_matrix = vcov(results))
# 多変量正規乱数を生成（β0, β1）
set.seed(123) # 再現性のため
samples = mvrnorm(n = 10000, mu = beta_hat, Sigma = vcov_matrix)
# 各サンプルから LD50 を計算
LD50_samples = -samples[,1] / samples[,2]
# 信頼区間（95%）
CI = quantile(LD50_samples, probs = c(0.025, 0.975))
# 結果表示
cat("シミュレーションによるLD50の95%信頼区間:\n")
print(CI)
862: 06/02(月)12:04 ID:GMuHFUYr(2/2) AAS
# データ
x <- c(-0.86, -0.3, -0.05, 0.73)
n <- c(5, 5, 5, 5)
y <- c(0, 1, 3, 5)
# JAGSモデル
model_string <- "
model {
for (i in 1:N) {
y[i] ~ dbin(p[i], n[i])
logit(p[i]) <- beta0 + beta1 * x[i]
}
# 事前分布（非情報的）
beta0 ~ dnorm(0.0, 0.001)
beta1 ~ dnorm(0.0, 0.001)
# LD50の定義
LD50 <- -beta0 / beta1
}
"
# JAGSに渡すデータ
data_jags <- list(
x = x,
n = n,
y = y,
N = length(y)
)
# 初期値
inits <- function() {
list(beta0 = rnorm(1, 0, 1), beta1 = rnorm(1, 0, 1))
}
# モデル作成と実行
model <- jags.model(textConnection(model_string), data = data_jags, inits = inits, n.chains = 3)
update(model, 1000) # バーンイン
# サンプリング
samples <- coda.samples(model, variable.names = c("beta0", "beta1", "LD50"), n.iter = 10000)
# 結果表示（LD50の95%信用区間）
summary(samples)
LD50_samples <- as.matrix(samples)[, "LD50"]
quantile(LD50_samples, probs = c(0.025, 0.975))
plot(samples)
863: 06/03(火)00:31 ID:CfA5PBxZ(1/2) AAS
x = c(-0.86, -0.3, -0.05, 0.73)
n = c(5, 5, 5, 5)
y = c(0, 1, 3, 5)
(data = cbind(x, n, y))
(response = cbind(y, n - y) )
results = glm(response ~ x, family = binomial(link="logit") )
summary(results)
b=MASS::mvrnorm(1e5,results$coef,S=vcov(results))
# log(p/(1-p)) = b1+b2*x
LD = function(b1,b2,p=0.5) (log(p/(1-p)) - b1)/b2
LD50=mapply(LD,b[,1],b[,2])
hist(LD50)
quantile(LD50,c(0.025,0.5,0.975))
HDInterval::hdi(LD50)
864: 06/03(火)23:40 ID:CfA5PBxZ(2/2) AAS
k=10
m=400
Nmax=1000
sim = \(){
N=sample(m:Nmax,1)
if(max(sample(1:N,k))==m) return(N)
}
N=NULL
while(length(N)<1000) N=c(N,sim())
summary(N)
quantile(N,c(0.5,0.95))
865: 06/04(水)07:49 ID:n75lIIio(1) AAS
set.seed(123)
solve = function(x,k=1e5){
f = function(x) runif(1,x-0.5,x+0.5)
y=replicate(k,mean(sapply(x,f)))
quantile(y,p=c(0.025,0.975))
}
solve(c(9,10,11,11,12))
set.seed(123) # 再現性のため
solve2 <- function(x, k = 1e5) {
# 各x_iに対して一様乱数を生成し、平均を計算（ブートストラップ）
bootstrap_means <- replicate(k, {
x_true <- runif(length(x), x - 0.5, x + 0.5)
mean(x_true)
})
# 95%信頼区間を計算
quantile(bootstrap_means, probs = c(0.025, 0.975))
}
# 実行例
x_rounded <- c(9, 10, 11, 11, 12)
solve2(x_rounded)
866: 06/05(木)13:25 ID:tGlaBVfa(1) AAS
> stancode(fit)
// generated with brms 2.22.0
functions {
/* compute monotonic effects
* Args:
* scale: a simplex parameter
* i: index to sum over the simplex
* Returns:
* a scalar between 0 and rows(scale)
*/
real mo(vector scale, int i) {
if (i == 0) {
return 0;
} else {
return rows(scale) * sum(scale[1:i]);
}
}
}
data {
int<lower=1> N; // total number of observations
array[N] int Y; // response variable
int<lower=1> K; // number of population-level effects
matrix[N, K] X; // population-level design matrix
int<lower=1> Kc; // number of population-level effects after centering
int<lower=1> Ksp; // number of special effects terms
int<lower=1> Imo; // number of monotonic variables
array[Imo] int<lower=1> Jmo; // length of simplexes
array[N] int Xmo_1; // monotonic variable
vector[Jmo[1]] con_simo_1; // prior concentration of monotonic simplex
int prior_only; // should the likelihood be ignored?
}
transformed data {
matrix[N, Kc] Xc; // centered version of X without an intercept
vector[Kc] means_X; // column means of X before centering
for (i in 2:K) {
means_X[i - 1] = mean(X[, i]);
Xc[, i - 1] = X[, i] - means_X[i - 1];
}
}
parameters {
vector[Kc] b; // regression coefficients
real Intercept; // temporary intercept for centered predictors
simplex[Jmo[1]] simo_1; // monotonic simplex
vector[Ksp] bsp; // special effects coefficients
}
transformed parameters {
real lprior = 0; // prior contributions to the log posterior
lprior += normal_lpdf(b[1] | 0.15, 0.3);
lprior += normal_lpdf(b[2] | 0.08, 0.3);
lprior += normal_lpdf(b[3] | 0.8, 0.3);
lprior += normal_lpdf(b[4] | 0.5, 0.3);
lprior += normal_lpdf(Intercept | -4, 2);
lprior += dirichlet_lpdf(simo_1 | con_simo_1);
lprior += normal_lpdf(bsp[1] | -0.5, 0.3);
}
model {
// likelihood including constants
if (!prior_only) {
// initialize linear predictor term
vector[N] mu = rep_vector(0.0, N);
mu += Intercept;
for (n in 1:N) {
// add more terms to the linear predictor
mu[n] += (bsp[1]) * mo(simo_1, Xmo_1[n]);
}
target += bernoulli_logit_glm_lpmf(Y | Xc, mu, b);
}
// priors including constants
target += lprior;
}
generated quantities {
// actual population-level intercept
real b_Intercept = Intercept - dot_product(means_X, b);
}
867: 06/06(金)04:41 ID:fR053ZqC(1/2) AAS
# ロジスティック回帰モデル
ACT = c(16,18,20,22,24,26,28)
n = c(2,7,14,26,13,14,3)
y = c(0,0,6,12,7,9,3)
fit = glm(cbind(y, n - y) ~ ACT, family = binomial())
#predict(fit, newdata = data.frame(ACT = 20)) |> plogis()
predict(fit, newdata = data.frame(ACT = 20)
,type="response")
cat("\n\n===== se.fit=TRUE =====\n\n")
pred1=predict(fit, newdata = data.frame(ACT = 20)
,type="response",se.fit=TRUE)
# 信頼区間（response))結果表示
ci=c(pred1$fit - 1.96*pred1$se.fit,pred1$fit + 1.96*pred1$se.fit)
cat("95% 信頼区間: [", ci[1], ",", ci[2], "]\n")
# 予測（log-odds とその標準誤差）
pred = predict(fit, newdata = data.frame(ACT = 20), se.fit = TRUE)
# 信頼区間（log-odds）
log_odds = pred$fit
se = pred$se.fit
lower_log_odds = log_odds - 1.96 * se
upper_log_odds = log_odds + 1.96 * se
# 確率（ロジスティック関数で変換）
prob = plogis(log_odds)
lower_prob = plogis(lower_log_odds)
upper_prob = plogis(upper_log_odds)
# 結果表示
cat("95% 信頼区間: [", lower_prob, ",", upper_prob, "]\n")
868: 06/06(金)04:43 ID:fR053ZqC(2/2) AAS
確率は0～1の範囲に制限されるため、直接計算は境界（0や1）に近い場合に不適切（例: 負の値や1超えの可能性）。また、ロジスティック回帰の標準誤差はログオッズスケールで計算されるため、確率スケールでの線形近似は精度が落ちる。
上下前次1-新書関写板覧索設栞歴
あと 134 ﾚｽあります
ｽﾚ情報 赤ﾚｽ抽出 画像ﾚｽ抽出 歴の未読ｽﾚ
ぬこの手 ぬこTOP 0.029s