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

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

765: １３２人目の素数さん [sage] 2025/03/25(火) 18:57:50.17 ID:JvJVqzfB

# 総列挙(３箱のうち２箱を表示）

ans1 # 1 1 8
m1=matrix(c("赤","黒","白","青")[ans1],ncol=2)
li1=cbind(m1[,1],'|',m1[,2]) |> noquote() ; li1

ans2 # 1 2 7
m2=matrix(c("赤","黒","白","青")[ans2],ncol=3)
li2=cbind(m2[,1],'|',m2[,2:3]) |> noquote()   ; li2

ans3 # 1 3 6
m3=matrix(c("赤","黒","白","青")[ans3],ncol=4)
li3=cbind(m3[,1],'|', m3[,2:4]) |> noquote()  ; li3

ans4 # 1 4 5
m4=matrix(c("赤","黒","白","青")[ans4],ncol=5)
li4=cbind(m4[,1],'|', m4[,2:5]) |> noquote()  ; li4

ans5 # 2 2 6
m5=matrix(c("赤","黒","白","青")[ans5],ncol=4)
li5=cbind(m5[,1:2],'|', m5[,3:4]) |> noquote() ; li5

ans6 # 2 3 5
m6=matrix(c("赤","黒","白","青")[ans6],ncol=5)
li6=cbind(m6[,1:2],'|', m6[,3:5]) |> noquote()  ; li6

ans7 # 2 4 4
m7=matrix(c("赤","黒","白","青")[ans7],ncol=8)
li7=cbind(m7[,1:4],'|', m7[,5:8]) |> noquote()  ; li7

ans8 # 3 3 4
m8=matrix(c("赤","黒","白","青")[ans8],ncol=6)
li8=cbind(m8[,1:3],'|', m8[,4:6]) |> noquote()  ; li8

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

766: １３２人目の素数さん [sage] 2025/03/25(火) 19:45:01.28 ID:JvJVqzfB

分け方の最初の30通りを列挙
[1,] 赤   |    赤   |    赤   赤   黒   黒   黒   白    白    青
[2,] 赤   |    黒   |    赤   赤   赤   黒   黒   白    白    青
[3,] 赤   |    白   |    赤   赤   赤   黒   黒   黒    白    青
[4,] 赤   |    青   |    赤   赤   赤   黒   黒   黒    白    白
[5,] 黒   |    黒   |    赤   赤   赤   赤   黒   白    白    青
[6,] 黒   |    白   |    赤   赤   赤   赤   黒   黒    白    青
[7,] 黒   |    青   |    赤   赤   赤   赤   黒   黒    白    白
[8,] 白   |    白   |    赤   赤   赤   赤   黒   黒    黒    青
[9,] 白   |    青   |    赤   赤   赤   赤   黒   黒    黒    白
[10,] 赤   |    赤   赤   |    赤   黒   黒   黒   白    白    青
[11,] 赤   |    赤   黒   |    赤   赤   黒   黒   白    白    青
[12,] 赤   |    赤   白   |    赤   赤   黒   黒   黒    白    青
[13,] 赤   |    赤   青   |    赤   赤   黒   黒   黒    白    白
[14,] 赤   |    黒   黒   |    赤   赤   赤   黒   白    白    青
[15,] 赤   |    黒   白   |    赤   赤   赤   黒   黒    白    青
[16,] 赤   |    黒   青   |    赤   赤   赤   黒   黒    白    白
[17,] 赤   |    白   白   |    赤   赤   赤   黒   黒    黒    青
[18,] 赤   |    白   青   |    赤   赤   赤   黒   黒    黒    白
[19,] 黒   |    赤   赤   |    赤   赤   黒   黒   白    白    青
[20,] 黒   |    赤   黒   |    赤   赤   赤   黒   白    白    青
[21,] 黒   |    赤   白   |    赤   赤   赤   黒   黒    白    青
[22,] 黒   |    赤   青   |    赤   赤   赤   黒   黒    白    白
[23,] 黒   |    黒   黒   |    赤   赤   赤   赤   白    白    青
[24,] 黒   |    黒   白   |    赤   赤   赤   赤   黒    白    青
[25,] 黒   |    黒   青   |    赤   赤   赤   赤   黒    白    白
[26,] 黒   |    白   白   |    赤   赤   赤   赤   黒    黒    青
[27,] 黒   |    白   青   |    赤   赤   赤   赤   黒    黒    白
[28,] 白   |    赤   赤   |    赤   赤   黒   黒   黒    白    青
[29,] 白   |    赤   黒   |    赤   赤   赤   黒   黒    白    青
[30,] 白   |    赤   白   |    赤   赤   赤   黒   黒    黒    青

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

767: １３２人目の素数さん [sage] 2025/03/25(火) 19:45:24.74 ID:JvJVqzfB

分け方の最後の30通りを列挙
[1,] 赤   黒   黒   |    赤   黒   青   |    赤   赤    白    白
[2,] 赤   赤   黒   |    黒   白   白   |    赤   赤    黒    青
[3,] 赤   赤   白   |    黒   黒   白   |    赤   赤    黒    青
[4,] 赤   黒   黒   |    赤   白   白   |    赤   赤    黒    青
[5,] 赤   黒   白   |    赤   黒   白   |    赤   赤    黒    青
[6,] 赤   赤   黒   |    黒   白   青   |    赤   赤    黒    白
[7,] 赤   赤   白   |    黒   黒   青   |    赤   赤    黒    白
[8,] 赤   赤   青   |    黒   黒   白   |    赤   赤    黒    白
[9,] 赤   黒   黒   |    赤   白   青   |    赤   赤    黒    白
[10,] 赤   黒   白   |    赤   黒   青   |    赤   赤    黒    白
[11,] 赤   赤   黒   |    白   白   青   |    赤   赤    黒    黒
[12,] 赤   赤   白   |    黒   白   青   |    赤   赤    黒    黒
[13,] 赤   赤   青   |    黒   白   白   |    赤   赤    黒    黒
[14,] 赤   黒   白   |    赤   白   青   |    赤   赤    黒    黒
[15,] 赤   黒   青   |    赤   白   白   |    赤   赤    黒    黒
[16,] 赤   黒   黒   |    黒   白   白   |    赤   赤    赤    青
[17,] 赤   黒   白   |    黒   黒   白   |    赤   赤    赤    青
[18,] 赤   白   白   |    黒   黒   黒   |    赤   赤    赤    青
[19,] 赤   黒   黒   |    黒   白   青   |    赤   赤    赤    白
[20,] 赤   黒   白   |    黒   黒   青   |    赤   赤    赤    白
[21,] 赤   黒   青   |    黒   黒   白   |    赤   赤    赤    白
[22,] 赤   白   青   |    黒   黒   黒   |    赤   赤    赤    白
[23,] 赤   黒   黒   |    白   白   青   |    赤   赤    赤    黒
[24,] 赤   黒   白   |    黒   白   青   |    赤   赤    赤    黒
[25,] 赤   黒   青   |    黒   白   白   |    赤   赤    赤    黒
[26,] 赤   白   白   |    黒   黒   青   |    赤   赤    赤    黒
[27,] 赤   白   青   |    黒   黒   白   |    赤   赤    赤    黒
[28,] 黒   黒   黒   |    白   白   青   |    赤   赤    赤    赤
[29,] 黒   黒   白   |    黒   白   青   |    赤   赤    赤    赤
[30,] 黒   黒   青   |    黒   白   白   |    赤   赤    赤    赤

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

768: １３２人目の素数さん [sage] 2025/03/25(火) 20:20:38.34 ID:JvJVqzfB

mat1=t(apply(ans1,1, \(x) c(x,fn_rest(v,x))))
M1=matrix(c("赤","黒","白","青")[mat1],ncol=10) |> noquote()
L1=cbind(M1[,1],'|',M1[,2],'|',M1[,3:10]) |> noquote() ; L1

mat2=t(apply(ans2,1, \(x) c(x,fn_rest(v,x))))
M2=matrix(c("赤","黒","白","青")[mat2],ncol=10) |> noquote()
L2=cbind(M2[,1],'|',M2[,2:3],'|',M2[,4:10]) |> noquote() ; L2

mat3=t(apply(ans3,1, \(x) c(x,fn_rest(v,x))))
M3=matrix(c("赤","黒","白","青")[mat3],ncol=10) |> noquote()
L3=cbind(M3[,1],'|',M3[,2:4],'|',M3[,5:10]) |> noquote() ; L3

mat4=t(apply(ans4,1, \(x) c(x,fn_rest(v,x))))
M4=matrix(c("赤","黒","白","青")[mat4],ncol=10) |> noquote()
L4=cbind(M4[,1],'|',M4[,2:5],'|',M4[,6:10]) |> noquote() ; L4

mat5=t(apply(ans5,1, \(x) c(x,fn_rest(v,x))))
M5=matrix(c("赤","黒","白","青")[mat5],ncol=10) |> noquote()
L5=cbind(M5[,1:2],'|',M5[,3:4],'|',M5[,5:10]) |> noquote() ; L5

mat6=t(apply(ans6,1, \(x) c(x,fn_rest(v,x))))
M6=matrix(c("赤","黒","白","青")[mat6],ncol=10) |> noquote()
L6=cbind(M6[,1:2],'|',M6[,3:5],'|',M6[,6:10]) |> noquote() ; L6

mat7=t(apply(ans7,1, \(x) c(x,fn_rest(v,x))))
M7=matrix(c("赤","黒","白","青")[mat7],ncol=10) |> noquote()
L7=cbind(M7[,1:2],'|',M7[,3:5],'|',M7[,6:10]) |> noquote() ; L7

mat8=t(apply(ans8,1, \(x) c(x,fn_rest(v,x))))
M8=matrix(c("赤","黒","白","青")[mat8],ncol=10) |> noquote()
L8=cbind(M8[,1:3],'|',M8[,4:6],'|',M8[,7:10]) |> noquote() ; L8

RES=rbind(L1,L2,L3,L4,L5,L6,L7,L8)|> noquote()
nrow(RES)
head(RES,30)
tail(RES,30)

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

771: １３２人目の素数さん [sage] 2025/03/31(月) 11:28:24.79 ID:u+Kd/O/2

# 赤玉a個、黒玉b個、白玉c個、青玉d個の合計(a+b+c+d)個の玉を空箱なしで3つの箱に分けて入れる。箱を区別しないとき、入れ方は何通りあるか？"

solve=function(a,b,c,d){
divide = function(n) {
indices = expand.grid(i = 0:n, j = 0:n)
indices = indices[indices$j >= indices$i, ]
result = lapply(1:nrow(indices), function(k) {
x = indices[k, ]
matrix(c(x$i, x$j - x$i, n - x$j), nrow = 1)
})
return(result)
}

reds = divide(a)
blacks = divide(b)
whites = divide(c)
blues = divide(d)
box3 = list()

combinations = expand.grid(red = reds, black = blacks, white = whites, blue = blues)

box3 = apply(combinations, 1, function(x) {
red = x[[1]]
black = x[[2]]
white = x[[3]]
blue = x[[4]]

box = do.call(rbind, lapply(1:3, function(i) {
sum_values = red[i] + black[i] + white[i] + blue[i]
if (sum_values > 0) c('red'=red[i], 'black'=black[i], 'white'=white[i], 'blue'=blue[i]) else NULL
}))

if (!is.null(box) & nrow(box) == 3) {
sorted_box = box[order(apply(box, 1, paste, collapse = ",")), ]
return(sorted_box)
}
return(NULL)
})

length(unique(Filter(Negate(is.null), box3)))
}

f=function(a,b,c,d){
ball=c(a,b,c,d)
x=(a+2)*(b+2)*(c+2)*(d+2)
y=(a+1)*(b+1)*(c+1)*(d+1)
z=floor(a/2+1)*floor(b/2+1)*floor(c/2+1)*floor(d/2+1)
u=ifelse(all(ball%%2==0),1,0)
v=ifelse(all(ball%%2==0),1,0)
x*y/96-(3*y-3*z-3*u+4*v)/6-u*v
}

p=sample(10,4)
p
solve(p[1],p[2],p[3],p[4])
f(p[1],p[2],p[3],p[4])

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

773: １３２人目の素数さん [sage] 2025/03/31(月) 11:35:19.58 ID:u+Kd/O/2

# 赤玉a個、黒玉b個、白玉c個、青玉d個の合計(a+b+c+d)個の玉を空箱なしで3つの箱に分けて入れる。箱を区別しないとき、入れ方は何通りあるか？"

reds = divide(a)
blacks = divide(b)
whites = divide(c)
blues = divide(d)
box3 = list()

combinations = expand.grid(red = reds, black = blacks, white = whites, blue = blues)

box3 = apply(combinations, 1, function(x) {
red = x[[1]]
black = x[[2]]
white = x[[3]]
blue = x[[4]]

if (!is.null(box) & nrow(box) == 3) {
sorted_box = box[order(apply(box, 1, paste, collapse = ",")), ]
return(sorted_box)
}
return(NULL)
})

length(unique(Filter(Negate(is.null), box3)))
}

p=sample(10,4)

solve(p[1],p[2],p[3],p[4])
f(p[1],p[2],p[3],p[4])

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

777: １３２人目の素数さん [sage] 2025/04/10(木) 06:36:28.17 ID:rUVTD6v0

rm(list=ls())
"
コイン9個を3行3列の枠に表を上にして配置する。
無作為に1個のコインを選ぶ。選ばれたコイン及び上下左右に隣接するコインを裏返す。
この操作をｎ回繰り返した後にすべてのコインが裏返っている確率をP[n]とする。
コインを選ぶ順番は区別する。同じコインを複数回選んでもよい。
"
#
Pn =\(n){
nu=(1-(-1)^n)*(6+4*3^n-4*5^n-7^n+9^n)
de = (2^9*9^n)
gmp::as.bigq(nu,de)
}
#

n2b <- function(num, N=2, digit = 9){
r=num%%N
q=num%/%N
while(q > 0 | digit > 1){
r=append(q%%N,r)
q=q%/%N
digit=digit-1
}
return(rev(r))
}
b2n=\(n) (n %*% 2^(0:8))[1,1]
M=matrix(0,ncol=512,nrow=512)
colnames(M)=0:511
rownames(M)=0:511
a1=c(1,1,0,1,0,0,0,0,0)
a2=c(1,1,1,0,1,0,0,0,0)
a3=c(0,1,1,0,0,1,0,0,0)
a4=c(1,0,0,1,1,0,1,0,0)
a5=c(0,1,0,1,1,1,0,1,0)
a6=c(0,0,1,0,1,1,0,0,1)
a7=c(0,0,0,1,0,0,1,1,0)
a8=c(0,0,0,0,1,0,1,1,1)
a9=c(0,0,0,0,0,1,0,1,1)
A=rbind(a1,a2,a3,a4,a5,a6,a7,a8,a9)
f1=\(i) apply(A,1,\(v) b2n(bitwXor(n2b(i),v)))
for(i in 0:511) M[as.character(i),as.character(f1(i))]=1/9
library(expm)
solve=\(n) (M %^%n)["0","511"]
solve=Vectorize(solve)
n=c(5,7,9,11,13,15,17,19,21)
cbind(n,'P[n]'=solve(n),Pn=as.numeric(Pn(n)))

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

778: １３２人目の素数さん [sage] 2025/04/10(木) 06:36:42.87 ID:rUVTD6v0

"
コイン9個を3行3列の枠に表を上にして配置する。コインを
１　２　３
４　?　６
７　８　９

とする。5のコインだけが裏返っており、残りは表が上である。
無作為に1個のコインを選ぶ。選ばれたコインはそのままで、上下左右に隣接するコインを裏返す。
同じコインを複数回選んでもよい。
この操作をｎ回繰り返した後にすべてのコインが裏返っている確率をQ[n]とする。
n=1,2,...,20の各々でQ[n]を求めよ。
"
B=A
for(i in 1:9) B[i,i]=0
f2=\(i) apply(B,1,\(v) b2n(bitwXor(n2b(i),v)))
M2=matrix(0,ncol=512,nrow=512)
colnames(M2)=0:511
rownames(M2)=0:511
for(i in 0:511) M2[as.character(i),as.character(f2(i))]=1/9
b2n(c(0,0,0,0,1,0,0,0,0)) # 16
solve2=\(n) (M2 %^%n)["16","511"]
solve2=Vectorize(solve2)
m=1:20
cbind(m,solve2(m))

#
flip2 <- function(v) {
# 隣接マップ
li <- list(
c(2, 4),        # 1
c(1, 3, 5),     # 2
c(2, 6),        # 3
c(1, 5, 7),   　# 4
c(2, 4, 6, 8),  # 5
c(3, 5, 9),     # 6
c(4, 8),        # 7
c(5, 7, 9),   　# 8
c(6, 8)         # 9
)

# 初期状態
state <- rep(0, 9)
state[5] <- 1

for (pos in v) {
# 選ばれた位置と隣接するコインを反転
for (idx in li[[pos]]) {
state[idx] <- 1 - state[idx]
}
}
return(state)
}
sim=\(n,verbose=FALSE){
v=sample(9,n,replace=TRUE)
res=sum(flip2(v))==9
if(verbose & res) cat(v,'\n')
invisible(res)
}
del=replicate(1e3,sim(3,TRUE))
res=numeric()
for(n in 1:20) res[n]=mean(replicate(1e5,sim(n)))
re=solve2(m)
cbind(m,round(re,4),res)
plot(re)
points(res,type='h',col=4,lwd=3)

Qn=solve2(1:100)
par(mar=c(4,4,4,4))
plot(Qn,type='l',xlab='n',ylab='Q[n]')
points(res)

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

779: １３２人目の素数さん [sage] 2025/04/10(木) 06:37:33.82 ID:rUVTD6v0

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

793: １３２人目の素数さん [sage] 2025/04/18(金) 12:53:56.52 ID:XpMQ9r6R

# 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)

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

794: １３２人目の素数さん [sage] 2025/04/23(水) 02:35:47.03 ID:t2PViPB2

# Candidate values for the number of red balls (0 to 100)
R_vals <- 0:100

# Observed data
k <- 4      # Number of red balls drawn
n <- 10     # Sample size
N <- 100    # Total number of balls

# Likelihood using the hypergeometric distribution
likelihood <- dhyper(k, R_vals, N - R_vals, n)

# Prior distribution: uniform
prior <- rep(1, length(R_vals))

# Unnormalized posterior
posterior_unnorm <- likelihood * prior

# Normalize to get the posterior distribution
posterior <- posterior_unnorm / sum(posterior_unnorm)

# MAP estimate (most probable value)
R_MAP <- R_vals[which.max(posterior)]

# Posterior mean (expected value)
R_mean <- sum(R_vals * posterior)

# 95% central credible interval
cumulative <- cumsum(posterior)
lower_CI <- R_vals[which(cumulative >= 0.025)[1]]
upper_CI <- R_vals[which(cumulative >= 0.975)[1]]

# 95% Highest Posterior Density Interval (HPDI)
sorted <- order(posterior, decreasing = TRUE)
cumsum_sorted <- cumsum(posterior[sorted])
HPDI_index <- sorted[which(cumsum_sorted <= 0.95)]
HPDI_range <- range(R_vals[HPDI_index])

# Display results
cat("MAP estimate:", R_MAP, "\n")
cat("Posterior mean:", round(R_mean, 2), "\n")
cat("95% central credible interval: [", lower_CI, ",", upper_CI, "]\n")
cat("95% HPDI: [", HPDI_range[1], ",", HPDI_range[2], "]\n")

# Plot posterior distribution
plot(R_vals, posterior, type = "h", lwd = 2,
main = "Posterior Distribution P(R | Data)",
xlab = "Number of Red Balls (R)", ylab = "Posterior Probability")
abline(v = c(lower_CI, upper_CI), col = "blue", lty = 2)
abline(v = HPDI_range, col = "red", lty = 3)
legend("topright", legend = c("95% Central CI", "95% HPDI"),
col = c("blue", "red"), lty = c(2,3))

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

795: １３２人目の素数さん [sage] 2025/04/24(木) 07:03:56.33 ID:9AuNSRyA

# 仮定
p_kokuritsu <- 0.01
p_f_ran <- 0.05
ratio_kokuritsu <- 0.1
ratio_f_ran <- 0.2
n_simulations <- 10000

# シミュレーション結果を格納するベクトル
kokuritsu_counts <- 0
f_ran_counts <- 0

for (i in 1:n_simulations) {
# ランダムに学歴を生成 (簡略化のため二択)
education <- sample(c("kokuritsu", "f_ran", "other"), 1, prob = c(ratio_kokuritsu, ratio_f_ran, 1 - ratio_kokuritsu - ratio_f_ran))

# 学歴に基づいて侮蔑語を使用するかどうかをシミュレート
uses_slur <- FALSE
if (education == "kokuritsu" && runif(1) < p_kokuritsu) {
uses_slur <- TRUE
kokuritsu_counts <- kokuritsu_counts + 1
} else if (education == "f_ran" && runif(1) < p_f_ran) {
uses_slur <- TRUE
f_ran_counts <- f_ran_counts + 1
}
}

# シミュレーション結果の表示
cat("シミュレーション回数:", n_simulations, "\n")
cat("難関国立大学卒で侮蔑語を使用した回数:", kokuritsu_counts, "\n")
cat("Fラン卒で侮蔑語を使用した回数:", f_ran_counts, "\n")

# 確率の比較 (あくまでシミュレーション上の数値)
prob_slur_kokuritsu <- kokuritsu_counts / (ratio_kokuritsu * n_simulations)
prob_slur_f_ran <- f_ran_counts / (ratio_f_ran * n_simulations)

cat("難関国立大学卒の人が侮蔑語を使う確率 (シミュレーション):", prob_slur_kokuritsu, "\n")
cat("Fラン卒の人が侮蔑語を使う確率 (シミュレーション):", prob_slur_f_ran, "\n")

if (prob_slur_f_ran > prob_slur_kokuritsu) {
cat("シミュレーションの結果では、Fラン卒の人の方が侮蔑語を使う可能性が高い傾向にあります。\n")
} else if (prob_slur_kokuritsu > prob_slur_f_ran) {
cat("シミュレーションの結果では、難関国立大学卒の人の方が侮蔑語を使う可能性が高い傾向にあります。\n")
} else {
cat("シミュレーションの結果では、両者の侮蔑語使用の可能性に大きな差は見られませんでした。\n")
}

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

800: １３２人目の素数さん [sage] 2025/04/29(火) 18:09:40.80 ID:pY4WJf3b

library(RcppAlgos)
library(fmsb)
library(matrixStats)

N <- 50

# Generate all combinations
cm <- comboGeneral(0:N, 3, repetition = FALSE)

# Pre-allocate
n_vec <- rep(N, 3)
success <- cm
failure <- matrix(n_vec, nrow = nrow(cm), ncol = 3, byrow = TRUE) - cm

# Modified fast Fisher function - correct implementation
fast_fisher <- function(success, failure) {
# Initialize p-value matrix
pvals <- matrix(NA_real_, nrow = nrow(success), ncol = 3)

# Perform pairwise comparisons
for (i in 1:nrow(success)) {
# 1 vs 2
pvals[i,1] <- fisher.test(matrix(c(success[i,1], failure[i,1],
success[i,2], failure[i,2]), nrow = 2))$p.value

# 1 vs 3
pvals[i,2] <- fisher.test(matrix(c(success[i,1], failure[i,1],
success[i,3], failure[i,3]), nrow = 2))$p.value

# 2 vs 3
pvals[i,3] <- fisher.test(matrix(c(success[i,2], failure[i,2],
success[i,3], failure[i,3]), nrow = 2))$p.value
}

# Bonferroni adjustment
pmin(pvals * 3, 1)  # Cap at 1 after adjustment
}

# Run with timing
system.time({
# Overall Fisher tests
overall_p <- apply(cbind(success, failure), 1, function(x) {
fisher.test(matrix(x, nrow = 2))$p.value
})

# Pairwise Fisher tests
pairwise_p <- fast_fisher(success, failure)
min_pairwise_p <- rowMins(pairwise_p, na.rm = TRUE)

# Filter condition
keep <- overall_p > 0.05 & min_pairwise_p < 0.05
result <- cm[keep, ]
})

# Print first few results
head(result)
nrow(result)  # Number of qualifying combinations

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

809: １３２人目の素数さん [sage] 2025/04/30(水) 08:07:44.58 ID:wedVH8wl

options(warn = -1)
alpha <- 0.05

sim_fisher <- function(N = 100) {
# Function to simulate data for three groups and perform Fisher's exact test.
# N: Total sample size.

# Determine sample sizes for each of the three groups.
A <- sample(1:(N - 2), 1) # Ensure A leaves enough room for B and C.
remaining <- N - A
if (remaining > 1) {
B <- sample(1:(remaining - 1), 1) # Ensure B is at least 1.
C <- N - A - B
ABC <- c(A, B, C) # Vector of group sizes.

# Randomly generate the number of successes for each group (cannot exceed group size).
abc <- sapply(ABC, function(x) if (x > 0) sample(0:x, 1) else 0)
x <- abc # Vector of number of successes per group.
n <- ABC # Vector of total samples per group.

# Create a data frame for the overall Fisher's exact test.
data_all <- data.frame(
group = factor(rep(c("A", "B", "C"), times = n)),
success = unlist(sapply(1:3, function(i) c(rep(1, x[i]), rep(0, n[i] - x[i]))))
)
table_all <- table(data_all$group, data_all$success)
fisher_pg <- fisher.test(table_all)$p.value # P-value of the overall Fisher's exact test.

# Perform pairwise Fisher's exact tests with Bonferroni correction.
pairwise_p_values <- numeric(3)
pairs <- combn(levels(data_all$group), 2, simplify = FALSE)

for (i in seq_along(pairs)) {
pair <- pairs[[i]]
subset_data <- subset(data_all, group %in% pair)
table_pair <- table(subset_data$group, subset_data$success)
pairwise_p_values[i] <- fisher.test(table_pair)$p.value # P-value of the pairwise Fisher's exact test.
}

min_pairwise_p_bonf <- min(p.adjust(pairwise_p_values, method = "bonferroni"), na.rm = TRUE) # Minimum Bonferroni-corrected p-value from pairwise tests.

list(fisher_pg = fisher_pg, min_pairwise_p_bonf = min_pairwise_p_bonf, x = x, n = n)
} else {
return(NULL) # Return NULL if group sizes are invalid.
}
}

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

810: １３２人目の素数さん [sage] 2025/04/30(水) 08:07:51.54 ID:wedVH8wl

# Find data where the overall Fisher's exact test is not significant,
# but at least one pairwise Fisher's exact test (with Bonferroni correction) is significant.
res_no_overall_sig_pairwise_sig <- NULL
while (is.null(res_no_overall_sig_pairwise_sig) || res_no_overall_sig_pairwise_sig$fisher_pg > alpha || res_no_overall_sig_pairwise_sig$min_pairwise_p_bonf > alpha) {
res_no_overall_sig_pairwise_sig <- sim_fisher()
}
cat("Data where overall Fisher's test is not significant, but pairwise is:\n")
print(res_no_overall_sig_pairwise_sig)
cat("\n")

# Find data where the overall Fisher's exact test is significant,
# but all pairwise Fisher's exact tests (with Bonferroni correction) are not significant.
res_overall_sig_no_pairwise_sig <- NULL
while (is.null(res_overall_sig_no_pairwise_sig) || res_overall_sig_no_pairwise_sig$fisher_pg < alpha || res_overall_sig_no_pairwise_sig$min_pairwise_p_bonf < alpha) {
res_overall_sig_no_pairwise_sig <- sim_fisher()
}
cat("Data where overall Fisher's test is significant, but pairwise is not:\n")
print(res_overall_sig_no_pairwise_sig)

options(warn = 0)

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

811: １３２人目の素数さん [sage] 2025/04/30(水) 08:12:32.33 ID:wedVH8wl

options(warn = -1)
alpha <- 0.05

sim_chisq <- function(N = 100) {
# Function to simulate data for three groups and perform Chi-squared test (without Yates' correction).
# N: Total sample size.

# Determine sample sizes for each of the three groups.
A <- sample(1:(N - 2), 1) # Randomly select a size for group A, ensuring space for B and C.
remaining <- N - A
if (remaining > 1) {
B <- sample(1:(remaining - 1), 1) # Randomly select a size for group B, ensuring space for C.
C <- N - A - B # Calculate the size for group C.
ABC <- c(A, B, C) # Vector containing the sample sizes of the three groups.

# Randomly generate the number of successes for each group (must be between 0 and the group size).
abc <- sapply(ABC, function(x) if (x > 0) sample(0:x, 1) else 0)
x <- abc # Vector containing the number of successes for each group.
n <- ABC # Vector containing the total number of trials for each group.

# Create a contingency table for the overall Chi-squared test.
contig_all <- rbind(s = x, f = n - x) # Rows: successes (s), failures (f); Columns: groups.
chisq_pg <- chisq.test(contig_all, correct = FALSE)$p.value # Perform Chi-squared test (no correction) and get the p-value.

# Perform pairwise proportion tests with Bonferroni correction.
pairwise_prop_p_values <- as.vector(
pairwise.prop.test(x, n, correct = FALSE, p.adj = "bon")$p.value
) # Perform pairwise proportion tests (no correction) and get Bonferroni-adjusted p-values.
min_pairwise_p_bonf <- min(pairwise_prop_p_values, na.rm = TRUE) # Get the minimum of the adjusted pairwise p-values.

# Return a list containing the overall p-value, the minimum Bonferroni-corrected pairwise p-value, successes, and total trials.
list(chisq_pg = chisq_pg, min_pairwise_p_bonf = min_pairwise_p_bonf, x = x, n = n)
} else {
return(NULL) # Return NULL if the group sizes are invalid.
}
}

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

812: １３２人目の素数さん [sage] 2025/04/30(水) 08:12:39.62 ID:wedVH8wl

# Find data where the overall Chi-squared test is not significant (p > alpha),
# but at least one pairwise proportion test (with Bonferroni correction) is significant (p <= alpha).
res_no_overall_sig_pairwise_sig <- NULL
while (is.null(res_no_overall_sig_pairwise_sig) || res_no_overall_sig_pairwise_sig$chisq_pg > alpha || res_no_overall_sig_pairwise_sig$min_pairwise_p_bonf > alpha) {
res_no_overall_sig_pairwise_sig <- sim_chisq() # Keep simulating until the condition is met.
}
cat("Data where overall Chi-squared test is not significant, but pairwise proportion test is:\n")
print(res_no_overall_sig_pairwise_sig)
cat("\n")

# Find data where the overall Chi-squared test is significant (p <= alpha),
# but all pairwise proportion tests (with Bonferroni correction) are not significant (p > alpha).
res_overall_sig_no_pairwise_sig <- NULL
while (is.null(res_overall_sig_no_pairwise_sig) || res_overall_sig_no_pairwise_sig$chisq_pg < alpha || res_overall_sig_no_pairwise_sig$min_pairwise_p_bonf < alpha) {
res_overall_sig_no_pairwise_sig <- sim_chisq() # Keep simulating until the condition is met.
}
cat("Data where overall Chi-squared test is significant, but pairwise proportion test is not:\n")
print(res_overall_sig_no_pairwise_sig)

options(warn = 0)

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

816: １３２人目の素数さん [sage] 2025/05/01(木) 09:09:32.28 ID:L1qIlz9/

Stanで作ったらコンパイルに時間がかかる。簡単なモデルはJAGSの方がいい。離散変数も扱えるし。
# JAGS model
library(rjags)

# Prepare the data
outcome_data <- c(rep(1, 17), rep(2, 21), rep(3, 15), rep(4, 21), rep(5, 20), rep(6, 6))
N <- length(outcome_data)
data_jags <- list(outcome = outcome_data, N = N)

# Initial values (adjust as needed)
inits_jags <- list(
list(alpha = rep(1, 6), eta = 1),
list(alpha = runif(6, 0.1, 2), eta = 5)
)

# Compile the model
model_jags <- jags.model(
file = "hierarchical_dice_model.jag",
data = data_jags,
n.chains = 2,
n.adapt = 1000
)

# Sampling
samples_jags <- coda.samples(
model = model_jags,
variable.names = c("prob_simplex", "alpha", "eta"),
n.iter = 4000
)

# Summary of the results
cat("\nJAGS Sampling Results Summary:\n")
summary(samples_jags)

# Extract posterior samples (prob_simplex)
prob_simplex_posterior_jags <- as.matrix(samples_jags[, grep("prob_simplex", varnames(samples_jags))])
head(prob_simplex_posterior_jags)

# Plotting (example: posterior distribution of probabilities for each outcome)
cat("\nPosterior Distribution Plots for Each Outcome:\n")
par(mfrow = c(2, 3))
for (i in 1:6) {
plot(prob_simplex_posterior_jags[, i], type = "l", main = paste("Prob[", i, "]"), xlab = "Iteration", ylab = "Probability")
abline(h = 1/6, col = "red", lty = 2)
}

dice_prob_mean=prob_simplex_posterior_jags
colors <- c("skyblue", "lightcoral", "lightgreen", "gold", "lightsalmon", "lightcyan")
for (i in 1:ncol(dice_prob_mean)) {
BEST::plotPost(dice_prob_mean[, i],
compVal=1/6,
xlab=paste("pip ", i),
xlim=c(0, 0.4),
main="",
col=colors[i],
border="black")
}

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

818: １３２人目の素数さん [sage] 2025/05/01(木) 11:18:32.55 ID:L1qIlz9/

ニッチな値の探索処理が終了しないコード

rm(list=ls())
library(fmsb)
library(parallel)

alpha <- 0.05

# Function to perform a single simulation
sim_single <- function(N = 100) {
A <- sample(1:N, 1)
while (A > N - 2) A <- sample(1:N, 1)
B <- sample(N - A - 1, 1)
C <- N - A - B
ABC <- c(A, B, C)
abc <- sapply(ABC, function(x) sample(x, 1))
x <- abc
n <- ABC
m <- rbind(s = x, f = n - x)

bonf_res <- pairwise.fisher.test(x, n, p.adj = 'bonf')
holm_res <- pairwise.fisher.test(x, n, p.adj = 'holm')
fdr_res <- pairwise.fisher.test(x, n, p.adj = 'fdr')
none_res <- pairwise.fisher.test(x, n, p.adj = 'none')

bonf <- min(bonf_res$p.value, na.rm = TRUE)
holm <- min(holm_res$p.value, na.rm = TRUE)
fdr <- min(fdr_res$p.value, na.rm = TRUE)
none <- min(none_res$p.value, na.rm = TRUE)

list(m = m, bonf = bonf, holm = holm, fdr = fdr, none = none)
}

# Function to find a result that meets the criteria using parallel processing
find_result_parallel_loop <- function(alpha = 0.05, num_cores = detectCores() - 1) {
# Create a cluster of worker processes
cl <- makeCluster(num_cores)
# Ensure the cluster is stopped when the function exits
on.exit(stopCluster(cl))

# Export the sim_single function to the worker processes
clusterExport(cl, "sim_single")

# Load the fmsb library in the worker processes
clusterEvalQ(cl, library(fmsb))

cat("Searching for a result that meets the criteria...\n")

while (TRUE) {
# Run simulations in parallel
results <- parLapply(cl, 1:num_cores, function(i) { # Run as many simulations as cores
sim_single()
})

# Check the results for the desired condition
for (res in results) {
if (res$bonf > alpha && res$holm < alpha) {
cat("Result meeting the criteria found:\n")
return(res)
}
}
# If no result found in this batch, continue the loop
}
}

# Find the result using parallel processing until found
res_parallel_loop <- find_result_parallel_loop(alpha = alpha)

# Output the result (will be printed within the loop when found)
print(res_parallel_loop)

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

822: １３２人目の素数さん [sage] 2025/05/02(金) 10:48:45.43 ID:+8QO9mMm

set.seed(123)
library(fmsb)

alpha <- 0.05
N <- 1000

sim <- function() {
# 群ごとのサンプルサイズ（不均等で可）
n1 <- sample(250:350, 1)  # 低用量
n2 <- sample(250:350, 1)  # 高用量
n3 <- N - n1 - n2         # プラセボ

# 各群の成功率設定（ペアでは差なし、プラセボ vs 合算で差あり）
p1 <- 0.30  # 低用量
p2 <- 0.31  # 高用量
p3 <- 0.22  # プラセボ（低め）

x1 <- rbinom(1, n1, p1)
x2 <- rbinom(1, n2, p2)
x3 <- rbinom(1, n3, p3)

m <- rbind(success = c(x1, x2, x3), failure = c(n1 - x1, n2 - x2, n3 - x3))

bonf_p <- suppressWarnings(min(as.vector(pairwise.fisher.test(m[1,], colSums(m), p.adj="bonf")$p.value), na.rm=TRUE))

combined <- matrix(c(x1 + x2, n1 + n2 - x1 - x2, x3, n3 - x3), nrow=2)
comb_p <- fisher.test(combined)$p.value

list(m = m, bonf = bonf_p, combined_p = comb_p, sizes = c(n1, n2, n3))
}

# 条件を満たすデータが出るまでループ
res <- sim()
while (!(res$bonf > alpha && res$combined_p < alpha)) {
res <- sim()
}

# 結果出力
res$m
cat("Bonferroni-adjusted pairwise p-value (min):", res$bonf, "\n")
cat("Combined treatment vs placebo p-value:", res$combined_p, "\n")
cat("Sample sizes:", res$sizes, "\n")

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

823: １３２人目の素数さん [sage] 2025/05/02(金) 11:02:18.75 ID:+8QO9mMm

set.seed(123)
library(fmsb)

alpha <- 0.05
N <- 1000

# Simulation function
sim <- function() {
# Random group sizes
n1 <- sample(250:350, 1)     # Low-dose
n2 <- sample(250:350, 1)     # High-dose
n3 <- N - n1 - n2            # Placebo
if (n3 < 100) return(NULL)   # Skip too-small placebo

# Success rates from uniform distribution
p1 <- runif(1)
p2 <- runif(1)
p3 <- runif(1)

# Binomial draws
x1 <- rbinom(1, n1, p1)
x2 <- rbinom(1, n2, p2)
x3 <- rbinom(1, n3, p3)

# 3-group matrix
m3 <- rbind(success = c(x1, x2, x3),
failure = c(n1 - x1, n2 - x2, n3 - x3))
colnames(m3) <- c("Low", "High", "Placebo")

# Add 4th group = combined (low + high)
x4 <- x1 + x2
n4 <- n1 + n2
m4 <- cbind(m3, Combined = c(x4, n4 - x4))

# Perform pairwise Fisher's exact tests across 4 groups
pw <- suppressWarnings(pairwise.fisher.test(m4[1,], colSums(m4), p.adj="bonf")$p.value)
pw_vals <- as.vector(pw)
pw_vals <- pw_vals[!is.na(pw_vals)]
names_all <- names(pw_vals)

# Identify significant pairs
sig_idx <- which(pw_vals < alpha)
sig_names <- names(pw_vals)[sig_idx]

# Check if only Combined vs Placebo is significant
is_valid <- length(sig_idx) == 1 &&
any(grepl("Placebo-Combined|Combined-Placebo", sig_names))

if (is_valid) {
return(list(m = m4, probs = c(p1, p2, p3), sizes = c(n1, n2, n3), pvals = pw))
} else {
return(NULL)
}
}

# Run until condition met
res <- NULL
while (is.null(res)) {
res <- sim()
}

# Output results
print(res$m)
cat("Success probabilities (Low, High, Placebo):", round(res$probs, 3), "\n")
cat("Sample sizes (Low, High, Placebo):", res$sizes, "\n")
cat("Pairwise Bonferroni-adjusted p-values:\n")
print(res$pvals)

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

824: １３２人目の素数さん [sage] 2025/05/02(金) 22:44:37.00 ID:056ygUN9

EMPAREG試験の解析をベイズでやっていたら、低用量高用量を統合する必要もなかったはず。
頻度主義統計で有意差がでない237:253の範囲でもプラセボよりイベント発生を抑制することが示せる。

############## 237:253 ################

# JAGSモデル文字列
model_string <- "
model {
for (i in 1:N) {
s[i] ~ dbin(theta[i], n[i])  # n[i] を別に与える
theta[i] ~ dbeta(1, 1)       # 非情報的事前分布
}
}
"

EMPA_REG=\(x,verbose=FALSE){
# データ
data_list <- list(
s = c(288, x, 490-x),
n = c(2333, 2345, 2342),
N = 3
)

# 初期値
# init_fun <- function() list(theta = runif(3, 0.05, 0.15))
init_fun <- function() list(theta = runif(3))

# モデル構築・初期化・実行
library(rjags)
model <- jags.model(textConnection(model_string),
data = data_list,
inits = init_fun,
n.chains = 3,
n.adapt = 1000)

# バーンインとサンプル取得
update(model, 1000)
samples <- coda.samples(model,
variable.names = c("theta"),
n.iter = 5000)
ms=as.matrix(samples)
Placebo=ms[,1]
Low=ms[,2]
High=ms[,3]
if(verbose){
ylim=c(0,max(max(density(Placebo)$y),  max(density(Low)$y),  max(density(High)$y)))
xu=max(max(density(Placebo)$x),  max(density(Low)$x),  max(density(High)$x))
xl=min(min(density(Placebo)$x),  min(density(Low)$x),  min(density(High)$x))
xlim=c(xl,xu)
plot(density(Placebo),xlim=xlim,ylim=ylim,xlab=quote(theta),ylab='',main='',col=8)
lines(density(Low),col='pink',lwd=2)
lines(density(High),col='red',lwd=2)
}
c(Low_Effective=mean(Low<Placebo),High_Effective=mean(High<Placebo) )
}

ans=sapply(237:253,EMPA_REG)
ans

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

825: １３２人目の素数さん [sage] 2025/05/04(日) 19:45:41.87 ID:Ie2Wyhjx

CRANからパッケージBESTが消えていたのでplotPostと同等機能の関数を復刻（不適当データ入力などエラー回避処理は面倒なのでやってない）。

https://i.imgur.com/GYBGnJs.png

This function helps visualize posterior distribution samples from Bayesian inference and displays various informative elements. It can show the mean, mode, median, a credible interval (either HDI or quantiles), the Region of Practical Equivalence (ROPE), and a comparison value.

Arguments:

posterior_samples: A numeric vector of posterior distribution samples.
credMass: The width of the credible interval (a numeric value between 0 and 1; default is 0.95).
ROPE: A numeric vector specifying the ROPE range (e.g., c(lower_bound, upper_bound)).
compVal: A numeric value for comparison.
showMean: A logical value indicating whether to display the mean (TRUE/FALSE).
showMode: A logical value indicating whether to display the mode (TRUE/FALSE).
showMedian: A logical value indicating whether to display the median (TRUE/FALSE).
showCurve: A logical value indicating whether to display a density plot (TRUE/FALSE; if FALSE, a histogram is shown).
showQuantile: A logical value indicating whether to display the credible interval as quantiles instead of HDI (TRUE/FALSE).
xlab: The label for the x-axis.
main: The title of the plot.
hist_color: The color of the histogram.
textSize: The size of the text elements.
yaxt: The y-axis style (default is to show only numbers).
...: Additional arguments to be passed to the plotting function.

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

826: １３２人目の素数さん [sage] 2025/05/05(月) 19:07:36.02 ID:LDY1RQtT

# This R function, riskratio.boot, calculates the risk ratio and its Highest Density Interval (HDI)
# using a bootstrap method. It takes the number of events and the total number of observations
# for two groups as input.

riskratio.boot <- function(r1, r2, n1, n2, nboot = 5000, conf.level = 0.95, verbose = FALSE){
# Combine the number of events and total observations for both groups.
r <- c(r1, r2)
n <- c(n1, n2)
# Create a 2x2 matrix representing the contingency table (number of events and non-events).
m <- cbind(r, n - r) ; m
# Perform bootstrap resampling to estimate the risk ratio distribution.
rr <- replicate(nboot, {
# Simulate the number of events in the first group by sampling with replacement
# from a population with the observed event rate.
sample(c(rep(1, r1), rep(0, n1 - r1)), n1, replace = TRUE) |> sum() -> R1
# Simulate the number of events in the second group similarly.
sample(c(rep(1, r2), rep(0, n2 - r2)), n2, replace = TRUE) |> sum() -> R2
# Calculate the risk ratio from the bootstrapped event counts.
(R1 / n1) / (R2 / n2)
})
# If verbose is TRUE, plot the posterior distribution of the risk ratio.
if(verbose){
source('plotpost.R') # Assuming 'plotpost.R' is a script for plotting posterior distributions.
plotpost(rr, compVal = 1, showCurve = 1, lwd = 4) # Plot the distribution with a comparison value of 1.
}
# Calculate the mean of the bootstrapped risk ratios.
b_mean <- mean(rr)
# Calculate the Highest Density Interval (HDI) of the bootstrapped risk ratios.
b_ci <- HDInterval::hdi(rr, credMass = conf.level)
# Return a list containing the bootstrap mean and the HDI.
list(b_mean = b_mean, b_ci = b_ci)
}

# Example usage of the riskratio.boot function.
riskratio.boot(244, 282, 2345, 2333)

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

827: １３２人目の素数さん [sage] 2025/05/05(月) 19:08:42.51 ID:LDY1RQtT

Description:

The `riskratio.boot` function in R estimates the risk ratio between two groups and provides its Highest Density Interval (HDI) using a bootstrap resampling approach. It takes the counts of events and the total number of observations for each of the two groups as input.

Usage:

riskratio.boot(r1, r2, n1, n2, nboot = 5000, conf.level = 0.95, verbose = FALSE)

Arguments:

* `r1`: The number of events in the first group.
* `r2`: The number of events in the second group.
* `n1`: The total number of observations in the first group.
* `n2`: The total number of observations in the second group.
* `nboot`: The number of bootstrap replicates to perform (default is 5000). A larger number of replicates generally provides more stable estimates.
* `conf.level`: The credibility level for the Highest Density Interval (HDI) (default is 0.95, corresponding to a 95% HDI).
* `verbose`: A logical value indicating whether to display a plot of the bootstrap distribution of the risk ratio (default is `FALSE`). If `TRUE`, it assumes a script named `plotpost.R` is available in the working directory for plotting.

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

828: １３２人目の素数さん [sage] 2025/05/05(月) 19:09:02.55 ID:LDY1RQtT

Details:

The function works by simulating the event outcomes in each group through bootstrap resampling. For each group, it draws `n1` (or `n2`) samples with replacement from a hypothetical population that has the observed proportion of events (`r1/n1` or `r2/n2`). The number of events in each resampled set (`R1` and `R2`) is then used to calculate a bootstrapped risk ratio `(R1/n1) / (R2/n2)`. This process is repeated `nboot` times to generate a distribution of risk ratios. The function then calculates the mean of this distribution and its Highest Density Interval (HDI), which represents the most credible range for the true risk ratio given the data and the bootstrap procedure.

If `verbose` is set to `TRUE`, the function attempts to plot the distribution of the bootstrapped risk ratios using a script named `plotpost.R`. This requires that the `plotpost.R` script exists in the current working directory and is capable of handling the vector of bootstrapped risk ratios.

Value:

The function returns a list with the following components:

* `b_mean`: The mean of the bootstrapped risk ratios.
* `b_ci`: The Highest Density Interval (HDI) of the bootstrapped risk ratios, as a vector with the lower and upper bounds.

Note:

The `verbose = TRUE` option depends on an external script `plotpost.R`. If you intend to use this option, ensure that the script is available and correctly implemented for plotting posterior-like distributions. The HDI is calculated using the `hdi` function from the `HDInterval` package, so this package must be installed (`install.packages("HDInterval")`) and loaded (`library(HDInterval)`) if you intend to use the default behavior.

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

831: １３２人目の素数さん [sage] 2025/05/10(土) 11:02:46.60 ID:ynDPH7B8

# ------------------------------------------------------------------------------
# ファイル名：logistic_regression_uraguchi_factors.R
# 目的：裏口入学の決定要因を評価するロジスティック回帰分析
# 考察対象の説明変数：学力、大学ランク (基準カテゴリ: A)、縁故、親の所得、寄付金
# ------------------------------------------------------------------------------

# データ生成 (大学ランクを因子型、基準レベル A)
set.seed(123)
n <- 1000
ranks_char <- sample(LETTERS[1:6], n, replace = TRUE)
ranks_factor <- factor(ranks_char, levels = LETTERS[1:6], ordered = TRUE) # 順序付きファクターとして生成

data <- data.frame(
裏口入学 = rbinom(n, 1, 0.2),
学力 = rnorm(n, mean = 50, sd = 10),
大学ランク = ranks_factor,
縁故 = rbinom(n, 1, 0.1),
親の所得 = rlnorm(n, meanlog = log(5000), sdlog = 0.3),
寄付金 = rlnorm(n, meanlog = log(100000), sdlog = 1.0)
)

# 支払額の生成
rank_numeric <- as.numeric(data$大学ランク) # A=1, B=2, ... 6
data$支払額 <- 10 * pmax(
300 + 100 * data$裏口入学 +
0.5 * (100 - data$学力) +
20 * rank_numeric +
50 * data$縁故 +
0.05 * data$親の所得 +
0.001 * data$寄付金 +
rnorm(n, mean = 0, sd = 50),
0
)

# ロジスティック回帰モデル (大学ランクが因子型として扱われ、基準カテゴリは A)
model_full <- glm(裏口入学 ~ 学力 + 大学ランク + 縁故 + 親の所得 + 寄付金, data = data, family = binomial)

# オッズ比と信頼区間の算出
confint_vals <- exp(confint(model_full))
odds_ratios_ci <- data.frame(
Variable = rownames(confint_vals)[-1],
CI_lower = confint_vals[-1, 1],
CI_upper = confint_vals[-1, 2]
)

# プロット用の変数ラベルを日本語化
label_map <- c(
"学力" = "学力",
"大学ランク.L" = "大学ランク B",
"大学ランク.Q" = "大学ランク C",
"大学ランク.C" = "大学ランク D",
"大学ランク^4" = "大学ランク E",
"大学ランク^5" = "大学ランク F",
"縁故" = "縁故",
"親の所得" = "親の所得",
"寄付金" = "寄付金"
)
odds_ratios_ci$日本語変数名 <- ifelse(odds_ratios_ci$Variable %in% names(label_map),
label_map[odds_ratios_ci$Variable],
odds_ratios_ci$Variable)

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

833: １３２人目の素数さん [sage] 2025/05/11(日) 21:13:09.76 ID:2CgV4g4d

# データ設定
n_placebo <- 1000; eff_placebo <- 24
n_old     <- 1000; eff_old     <- 40
n_new     <- 1000; eff_new     <- 25

# 有効率
p_placebo <- eff_placebo / n_placebo
p_old     <- eff_old     / n_old
p_new     <- eff_new     / n_new

# 比較：旧薬 vs 偽薬（有意差）
m1 <- matrix(c(eff_old, n_old - eff_old, eff_placebo, n_placebo - eff_placebo), nrow = 2)
test1 <- prop.test(m1, correct = FALSE)

# 比較：新薬 vs 偽薬（有意差なし）
m2 <- matrix(c(eff_new, n_new - eff_new, eff_placebo, n_placebo - eff_placebo), nrow = 2)
test2 <- prop.test(m2, correct = FALSE)

# 比較：旧薬 vs 新薬（非劣性検定）
# 非劣性マージン
M <- -0.10
# 差（新薬 - 旧薬）
diff <- p_new - p_old
# 標準誤差（差の95%信頼区間に使用）
se <- sqrt(p_new*(1 - p_new)/n_new + p_old*(1 - p_old)/n_old)
z <- qnorm(0.025, lower.tail = FALSE)
lower_CI <- diff - z * se

# 非劣性判定
non_inferior <- lower_CI > M

# 結果表示
cat("=== 旧薬 vs 偽薬 ===\n")
print(test1)
cat("\n=== 新薬 vs 偽薬 ===\n")
print(test2)
cat("\n=== 非劣性検定（旧薬 vs 新薬） ===\n")
cat(sprintf("差（新薬 - 旧薬） = %.3f\n", diff))
cat(sprintf("95%% CI = [%.3f, %.3f]\n", diff - z*se, diff + z*se))
cat(sprintf("非劣性マージン = %.3f\n", M))
cat(sprintf("非劣性判定: %s\n", ifelse(non_inferior, "非劣性あり", "非劣性なし")))

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

834: １３２人目の素数さん [sage] 2025/05/13(火) 13:21:07.90 ID:L+Wotuil

ド底辺シリツ医大の三法則を与えたらAIが１２法則まで拡張してくれました。

ド底辺医大の十二箇条 (The Twelve Laws of Do-Teihen Medical School)

第1法則
ド底辺シリツ医大が悪いのではない、本人の頭が悪いんだ。
It is not the bottom medical school but its enrollee that is despicable, which deserves to be called a bona fide moron beyond redemption.

第2法則
ド底辺シリツ医大卒は恥ずかしくて、学校名を皆さま言いません。
The graduates of Do-Teihen are so ashamed that none of them dare to mention their own alma mater they had gone through.

第3法則
ド底辺特殊シリツ医大卒は裏口入学の負い目から裏口馬鹿を暴く人間を偽医者扱いしたがる。
The Do-Teihen graduates are so ashamed of having bought their way into the exclusively bottom-leveled medical school that they tend to call a genuine doctor a charlatan who elucidates their imbecility.

第4法則
ド底辺医大卒は、偏差値や出身校を気にするなと言いながら、自分の子どもには絶対にそんな大学へは行かせたくないと思っている。
While claiming that academic ranking or alma mater does not matter, a Do-Teihen graduate would never let their own child attend such a university.

第5法則
ド底辺医大卒は、裏口入学を否定しない。否定できない。なぜなら、実際に自分がそうだからである。
A Do-Teihen graduate never denies the existence of backdoor admissions—because deep down, they know they were one of them.

第6法則
ド底辺医大は、「差別するな」と叫びながら、偏差値・財力・コネがない者を最もあからさまに差別する。
While crying out against discrimination, Do-Teihen medical schools are the very institutions that blatantly discriminate against those without test scores, wealth, or connections.

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

835: １３２人目の素数さん [sage] 2025/05/13(火) 13:21:14.83 ID:L+Wotuil

第7法則
ド底辺医大卒は「実力で入った」と言うが、その“実力”の定義を決して口にしない。
A Do-Teihen graduate may claim, “I got in on merit,” but they will never define what that 'merit' actually was.

第8法則
ド底辺医大卒の最大の敵は、同級生ではなく、偏差値という現実である。
The greatest enemy of a Do-Teihen graduate is not their classmates—but the cold, numerical reality of standardized test scores.

第9法則
ド底辺医大では、人格者は浮く。媚びる者とカネ持ちが残る。
In Do-Teihen med schools, the virtuous are outcasts; only flatterers and the wealthy thrive.

第10法則
ド底辺医大に入る者は、嘘をついて入学し、嘘をつき続けて卒業する。
Those who enter Do-Teihen medical schools do so with lies—and graduate by continuing to lie.

第11法則
ド底辺医大卒は、知性を持つ批判者を最も憎む。それは、自分が決してなれない姿だからだ。
Graduates of Do-Teihen reserve their deepest hatred for intelligent critics—because those critics reflect everything they can never become.

第12法則
ド底辺医大は、医者を育てる場ではない。医師免許を与える「通行証発行所」である。
Do-Teihen is not a school to train doctors. It is a toll booth that issues medical licenses for a price.

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

837: １３２人目の素数さん [sage] 2025/05/16(金) 16:37:03.25 ID:s89ybxV8

イベント発生が人数比で
臨床試験１で　旧薬 vs プラセボで　5/201 vs 19/202
臨床試験２で　新薬 vs　旧薬 で　  9/203 vs  5/204
であったとき
(1) 新薬がプラセボより劣る確率を計算せよ。
(2) 新薬はプラセボより有意差をもって有効といえるか？
計算に必要な条件は適宜設定してよい。
例：イベント発生は独立事象である

library(rjags)
library(coda)

worth_than_placebo <- function(r0, n0, r1, n1, r2, n2, r3, n3){
model_string <- '
model {
# 試験 (旧薬 vs 偽薬)
r1 ~ dbin(p1, n1)
p1 ~ dbeta(1, 1)
r0 ~ dbin(p0, n0)
p0 ~ dbeta(1, 1)

# 試験 (新薬 vs 旧薬)
r2 ~ dbin(p2, n2)
p2 ~ dbeta(1, 1)
r3 ~ dbin(p3, n3)
p3 ~ dbeta(1, 1)

# parameters
p2_est <- p2
p0_est <- p0
p2_worse_than_p0 <- step(p2_est - p0_est)
}
'
data <- list(r1=r1 , n1=n1 , r0=r0 ,n0=n0, r2=r2 , n2=n2 , r3=r3 , n3=n3)
jags_model <- jags.model(file=textConnection(model_string), data=data, n.chains=3,
n.adapt=3000, quiet = TRUE)
update(jags_model, n.iter=2000)
jags_samples <- coda.samples(jags_model, variable.names=c("p2_est", "p0_est", "p2_worse_than_p0"),
n.iter=10000, thin=1)
summary(jags_samples)
js <- as.data.frame(as.matrix(jags_samples))
names(js)

source("plotpost.R")
layout(matrix(c(1,2,3,3), 2, 2, byrow=TRUE))
plotpost(js$p2_est, col='lightcoral',xlab="新薬",cex.lab=1.5,main="")
plotpost(js$p0_est, col='lightgray', xlab="プラセボ",cex.lab=1.5,main="")
plotpost(js$p0_est - js$p2_est, compVal = 0, col=c('lightcoral', 'lightgray'),
xlab="プラセボ - 新薬", main="",cex.lab=1.5)
HDInterval::hdi(js$p0_est - js$p2_est) |> print()
mean(js$p2_worse_than_p0)
}

result <- worth_than_placebo(r0=19, n0=202, r1=5,  n1=201, r2=9,  n2=203, r3=5,  n3=204)
print(paste("新薬がプラセボより劣る確率:", result))

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

838: １３２人目の素数さん [sage] 2025/05/17(土) 02:21:13.42 ID:zAzyVzie

>>837
# --- 必要パッケージ ---
library(rjags)
library(coda)
library(HDInterval)

# --- データ定義 ---
data_list <- list(
r0 = 19, n0 = 202,  # プラセボ
r1 = 5,  n1 = 201,  # 旧薬（試験1）
r2 = 9,  n2 = 203,  # 新薬
r3 = 5,  n3 = 204   # 旧薬（試験2）
)

# --- 階層モデル定義 ---
model_hier <- "
model {
r0 ~ dbin(p0, n0)
r1 ~ dbin(p1, n1)
r2 ~ dbin(p2, n2)
r3 ~ dbin(p3, n3)

p0 ~ dbeta(1, 1)
p2 ~ dbeta(1, 1)

mu_old ~ dbeta(1, 1)
tau ~ dgamma(0.001, 0.001)  # 弱情報事前分布
p1 ~ dbeta(mu_old * tau, (1 - mu_old) * tau)
p3 ~ dbeta(mu_old * tau, (1 - mu_old) * tau)

p2_worse_than_p0 <- step(p2 - p0)

rd_p0_p1 <- p0 - p1
rd_p1_p2 <- p1 - p2
rd_p0_p3 <- p0 - p3

}
"

jags_model <- jags.model(textConnection(model_hier),
data = data_list, n.chains = 2, quiet=TRUE)
update(jags_model, 3000, progress.bar="none")
jags_samples <- coda.samples(jags_model,
c("p0","p1","p2","p3",
"p2_worse_than_p0", "rd_p0_p1","rd_p1_p2", "rd_p0_p3"),
n.iter=10000, progress.bar="none")
gelman.plot(jags_samples)
plot(jags_samples)

js <- as.data.frame(as.matrix(jags_samples))
mean(js$p2_worse_than_p0)
hdi(js$rd_p0_p1) # 旧薬(試験1)  vs プラセボ
hdi(js$rd_p1_p2) # 旧薬(試験2)  vs 新薬
hdi(js$rd_p0_p3) # 旧薬(試験2)  vs プラセボ　仮想
hdi(js$p0-js$p2) # 新薬　　　　 vs プラセボ　仮想

source("plotpost.R")
layout(matrix(c(1,2,3,3), 2, 2, byrow=TRUE))
plotpost(js$p2, col='lightcoral',xlab="新薬",cex.lab=1.5,main="")
plotpost(js$p0, col='lightgreen', xlab="プラセボ",cex.lab=1.5,main="")
plotpost(js$p0 - js$p2, compVal = 0, col=c('lightcoral', 'lightgreen'),
xlab="プラセボ - 新薬", cex.main=2,main="二項分布階層モデル",cex.lab=1.5, breaks="scott")

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

839: １３２人目の素数さん [sage] 2025/05/17(土) 07:49:04.94 ID:Vpav5/5q

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

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

840: １３２人目の素数さん [sage] 2025/05/20(火) 23:31:13.05 ID:gwaBTE4C

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)

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

842: １３２人目の素数さん [sage] 2025/05/24(土) 02:35:14.53 ID:VetM3rz7

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

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

843: １３２人目の素数さん [sage] 2025/05/24(土) 08:17:52.86 ID:VetM3rz7

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

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

844: １３２人目の素数さん [sage] 2025/05/24(土) 08:18:19.92 ID:VetM3rz7

# 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")

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

845: １３２人目の素数さん [sage] 2025/05/24(土) 08:37:05.75 ID:VetM3rz7