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

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

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-新書関写板覧索設栞歴

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

ぬこの手ぬこTOP 0.011s