高校数学の質問スレ（医者・東大卒専用） Part438

高校数学の質問スレ（医者・東大卒専用） Part438 (991ﾚｽ)
上下前次 1-新
通常表示 512ﾊﾞｲﾄ分割ﾚｽ栞

857: １３２人目の素数さん [sage] 2025/05/31(土) 05:25:02.62 ID:jzcOJBMt

#' @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)
)
}

http://rio2016.5ch.net/test/read.cgi/math/1723152147/857

866: １３２人目の素数さん [sage] 2025/06/05(木) 13:25:27.68 ID:tGlaBVfa

> 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);
}

http://rio2016.5ch.net/test/read.cgi/math/1723152147/866

867: １３２人目の素数さん [sage] 2025/06/06(金) 04:41:04.79 ID:fR053ZqC

# ロジスティック回帰モデル
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")

http://rio2016.5ch.net/test/read.cgi/math/1723152147/867

871: １３２人目の素数さん [sage] 2025/06/12(木) 09:46:29.84 ID:z7P0Lqdi

Bayesian Computation with RでRのコードが理解できなかった。バグだとおもったのだが、
https://bayesball.github.io/bcwr/corrections.2nd.edition.txt
のerrataにも掲載がないのでAIに聞いてみた。

>>
対数ヤコビアン項が間違っていると思う。

# theta=c(log(eta/(1-eta)),log(K))
> LearnBayes::betabinexch
function (theta, data)
{
eta = exp(theta[1])/(1 + exp(theta[1]))
K = exp(theta[2])
y = data[, 1]
n = data[, 2]
N = length(y)
logf = function(y, n, K, eta) lbeta(K * eta + y, K * (1 -
eta) + n - y) - lbeta(K * eta, K * (1 - eta))
val = sum(logf(y, n, K, eta))
val = val + theta[2] - 2 * log(1 + exp(theta[2]))
return(val)
}
<bytecode: 0x000001a5a980e758>
<environment: namespace:LearnBayes>

これが正しいのでは？
betabinexch <- function (theta, data)
{
eta = exp(theta[1])/(1 + exp(theta[1]))
K = exp(theta[2])
y = data[, 1]
n = data[, 2]
N = length(y)
logf = function(y, n, K, eta)
lbeta(K * eta + y, K * (1 - eta) + n - y) - lbeta(K * eta, K * (1 - eta))
val = sum(logf(y, n, K, eta))
val = val + theta[1] + theta[2] - 2 * log(1 + exp(theta[1])) # log Jacobian term
return(val)
}
<<

いずれのAIも
>あなたの指摘は正しいです。対数ヤコビアン項に問題があります。
という趣旨の返事が返ってきた。

http://rio2016.5ch.net/test/read.cgi/math/1723152147/871

879: １３２人目の素数さん [sage] 2025/06/15(日) 08:32:56.90 ID:MIIBNstg

pdf2hdi <- function(pdf,
xMIN=0,
xMAX=1,
cred=0.95,
Print=TRUE,
nxx=1001){

xx=seq(xMIN,xMAX,length=nxx)
xx=xx[-nxx]
xx=xx[-1]
xmin=xx[1]
xmax=xx[nxx-2]

AUC=integrate(pdf,xmin,xmax)$value
PDF=function(x)pdf(x)/AUC

cdf <- function(x) integrate(PDF,xmin,x)$value

ICDF <- function(x) uniroot(function(y) cdf(y)-x,c(xmin,xmax))$root
ICDF=Vectorize(ICDF)

hdi=HDInterval::hdi(ICDF,credMass=cred)

print(c(hdi[1],hdi[2]),digits=5)

if(Print){
par(mfrow=c(3,1))

plot(xx,sapply(xx,PDF),main='pdf',type='h',xlab='x',ylab='Density',col='lightgreen',bty='l')
legend('top',bty='n',legend=paste('HDI:',round(hdi,3)))

plot(xx,sapply(xx,cdf),main='cdf',type='h',xlab='x',ylab='Probability',col='lightblue',bty='l')

pp=seq(0,1,length=nxx)
pp=pp[-nxx]
pp=pp[-1]
plot(pp,sapply(pp,ICDF),type='l',xlab='p',ylab='x',main='ICDF',bty='l')

par(mfrow=c(1,1))
}

invisible(ICDF)
}
ICDF=pdf2hdi(function(x) dbeta(x,2,5))
hist(ICDF(seq(1e-12,1-1e-12,le=1000)))

AIの評価
まとめ
この pdf2hdi 関数は、数値積分と数値最適化 (uniroot) を巧みに組み合わせることで、任意のPDFからICDFを頑健に導出し、さらに統計的な要約であるHDIを計算する、非常に実用的かつ教育的なコードです。両端での数値計算の回避や正規化といった細部への配慮も素晴らしいです。
これにより、複雑な確率分布でも、そこからサンプリングしたり、HDIを求めたりといった解析が可能になります。

http://rio2016.5ch.net/test/read.cgi/math/1723152147/879

914: １３２人目の素数さん [sage] 2025/06/22(日) 11:18:18.20 ID:x1q6L9ud

Rだと小数に表示なのでMathematicaで計算。

80353937215217784622318561352545314301219984594089748478819
Out[4]= {------------------------------------------------------------,
803469022129495137770981046170581301261101496891396417650688

20088148412335002790621607242593839987034997410025084906975
>    ------------------------------------------------------------,
100433627766186892221372630771322662657637687111424552206336

40174538074102937850551109060704760584515594716727372518775
>    ------------------------------------------------------------,
200867255532373784442745261542645325315275374222849104412672

20086182069423287937895596774418478938984256623118935282825
>    ------------------------------------------------------------,
100433627766186892221372630771322662657637687111424552206336

160682421553118032860967525489926293991268384984698970271425
>    ------------------------------------------------------------,
803469022129495137770981046170581301261101496891396417650688

1255310425166894422771138921081126550349948690321031989171
>    -----------------------------------------------------------}
12554203470773361527671578846415332832204710888928069025792

In[5]:= N@ans1

Out[5]= {0.100009, 0.200014, 0.200005, 0.199995, 0.199986, 0.0999912}

In[6]:= Flatten@Position[ans1,Max[ans1]]

Out[6]= {2}

In[7]:=

http://rio2016.5ch.net/test/read.cgi/math/1723152147/914

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

あと 54 ﾚｽあります
ｽﾚ情報赤ﾚｽ抽出画像ﾚｽ抽出歴の未読ｽﾚ AAｻﾑﾈｲﾙ

ぬこの手ぬこTOP 0.011s