高校数学の質問スレ（医者・東大卒専用） Part438 (991ﾚｽ)
上下前次1-新
抽出解除 必死ﾁｪｯｶｰ(本家) (べ) 自ID ﾚｽ栞 あぼーん

---

847(1): １３２人目の素数さん [sage] 05/25(日)04:22 ID:P4nhnL8B(1/6)
# 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: １３２人目の素数さん [sage] 05/25(日)05:57 ID:P4nhnL8B(2/6)
＞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: １３２人目の素数さん [sage] 05/25(日)06:42 ID:P4nhnL8B(3/6)
>>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: １３２人目の素数さん [sage] 05/25(日)06:42 ID:P4nhnL8B(4/6)
# 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: １３２人目の素数さん [sage] 05/25(日)06:51 ID:P4nhnL8B(5/6)
今回スムーズに機能拡張できた理由は、以下のようにプロトタイプが非常に明快かつ健全だったことが大きな要因です。

✅ プロトタイプの良さが拡張性を支えた理由
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: １３２人目の素数さん [sage] 05/25(日)10:38 ID:P4nhnL8B(6/6)
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)
}

---

上下前次1-新書関写板覧索設栞歴

---

ｽﾚ情報 赤ﾚｽ抽出 画像ﾚｽ抽出 歴の未読ｽﾚ AAｻﾑﾈｲﾙ

---

ぬこの手 ぬこTOP 1.695s*