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

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

このｽﾚｯﾄﾞは過去ﾛｸﾞ倉庫に格納されています｡
次ｽﾚ検索歴削→次ｽﾚ栞削→次ｽﾚ過去ﾛｸﾞﾒﾆｭｰ

753: １３２人目の素数さん [sage] 2025/03/12(水) 21:31:17.69 ID:pDKAHj/P

「鹿児島 私立中学 寮 事件」や「学校名 寮 いじめ」で画像検索したらどこかわかるけど、絶対に探すなよ
i.imgur.com/s0F0tJb.jpeg i.imgur.com/4nmblMW.png

ソース 鹿児島放送 動画 https://www.tiktok.com/@kkb_kagoshima5/video/7428406277253926162

【鹿児島】私立中男子が入浴中に寮で大けが…「加害生徒の処分を」被害生徒の父が会見、いじめ重大事態に認定した学校へ調査と適切な対応求める

鹿児島県内の私立中学校の寮で、入浴していた１年生の男子生徒が、同学年の複数の男子生徒に両足首を引っ張られて転倒し、右手にまひが残るけがを負っていたことが７日、分かった。

学校はいじめ防止対策推進法に基づく重大事態に認定。今後、第三者委員会を設置して調査を進める。

被害生徒の保護者や学校によると、発生は２０２４年６月２５日。寮の共同浴場で、一緒に入浴していた男子生徒に両足首を引っ張り上げられた。転倒した際、湯を張った浴槽の床に後頭部と首を強く打ち付けた。

７月上旬、別の男子生徒に拳で肩をたたかれて痛みが強くなったため、１０日に養護教諭に相談し発覚。同日に病院を受診し、頸髄不全損傷（全治約１カ月）と診断された。

保護者は８月中旬、県警に被害届を提出し、受理されている。

９月末、いじめ認定を求める申立書が被害生徒の保護者から出され、対応を再協議。１０月３日に重大事態に認定した。第三者委員会の設置時期は未定。
教頭は「重く受け止めている。症状が早く改善し、普段通り学校生活を送れるように対応する」としている。

この問題を受け、４０代の父親が２１日、鹿児島市の県庁で会見を開いた。「学校は、加害生徒への処分を検討してほしい」と訴えた。父親は「大けがなのに相応の処分がないのは納得できない。しっかりと調査して、適切に対応してほしい」と求めた。

南日本新聞
373news.com/_news/storyid/202380/
373news.com/_news/storyid/203141/

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

754: １３２人目の素数さん [sage] 2025/03/19(水) 16:01:22.42 ID:TUD6PvYL

"
6人でババ抜きをした時の初期状態の手札の枚数の期待値は?
カードは52枚+ジョーカー
"

sim=function(){
# 1から13までの数字が各々4枚と0を含む53枚のリストを作成
numbers <- rep(c(0, 1:13), times = c(1, rep(4, 13)))
# 数字をシャッフルする
shuffled_numbers <- sample(numbers)
# 9枚を5組、8枚を1組に分ける
group_size_9 <- rep(9, 5)  # 9枚のグループ5つ
group_size_8 <- 8  # 8枚のグループ1つ
group_sizes <- c(rep(9, 5), 8)  # 9枚5組 + 8枚1組
# グループに分ける
groups <- split(shuffled_numbers, rep(1:6, times = group_sizes))
# ステップ4: 各グループで同じ数字があれば2個を組み合わせて捨てる
discarded_numbers <- list()  # 捨てられる数字を保存
remaining_numbers <- list()  # 残る数字を保存
for (i in 1:length(groups)) {
group <- groups[[i]]
count <- table(group)  # 各数字の出現回数をカウント
discarded <- count %/% 2  # 同じ数字のペアを捨てる（整数除算）
remaining <- count %% 2    # 残りの数字（1枚残るもの）
# 捨てた数字と残りの数字をリストに保存
discarded_numbers[[i]] <- rep(names(discarded)[discarded > 0], discarded[discarded > 0])
remaining_numbers[[i]] <- rep(names(remaining)[remaining > 0], remaining[remaining > 0])
}
# 残りの枚数を計算
length(unlist(remaining_numbers))
}

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

756: １３２人目の素数さん [sage] 2025/03/24(月) 08:15:20.18 ID:Wji9obyD

赤玉5個、黒玉3個、白玉2個を(4個, 3個, 3個)に分ける方法（もらう人を区別しない）」を全て列挙せよ

library(gtools)  # permutations を使用

# 各色の玉の個数
balls <- c(赤 = 4, 黒 = 3, 白 = 2, 青 = 1)

# 各箱のサイズ
box_sizes <- c(4, 3, 3)

# 指定したサイズで玉を選ぶ関数
generate_combinations <- function(balls, size) {
colors <- names(balls)
all_combinations <- permutations(n = length(colors), r = size, v = colors, repeats.allowed = TRUE)

# 各組み合わせが元の玉の数を超えないようにフィルタリング
valid_combinations <- all_combinations[apply(all_combinations, 1, function(combo) {
all(table(combo) <= balls[names(table(combo))])
}), ]

return(unique(t(apply(valid_combinations, 1, sort))))  # 各組み合わせをソートして重複を削除
}

# すべての分割を求める
unique_distributions <- list()

first_box_choices <- generate_combinations(balls, box_sizes[1])

for (first_box in 1:nrow(first_box_choices)) {
remaining1 <- balls
first_selection <- first_box_choices[first_box, ]

for (color in first_selection) {
remaining1[color] <- remaining1[color] - 1
}

second_box_choices <- generate_combinations(remaining1, box_sizes[2])

for (second_box in 1:nrow(second_box_choices)) {
remaining2 <- remaining1
second_selection <- second_box_choices[second_box, ]

for (color in second_selection) {
remaining2[color] <- remaining2[color] - 1
}

# 第3箱は残りの玉で確定
third_box <- rep(names(remaining2), times = remaining2)

# 各箱内をソートし、3つの箱の並びも統一
sorted_distribution <- list(sort(first_selection), sort(second_selection), sort(third_box))
sorted_distribution <- sorted_distribution[order(sapply(sorted_distribution, paste, collapse = ""))]  # 箱順も統一

# **文字列化してリストを扱いやすくする**
unique_distributions <- append(unique_distributions, list(paste(sapply(sorted_distribution, paste, collapse = ","), collapse = " | ")))
}
}

# 重複を削除
unique_distributions <- unique(unique_distributions)

# 総数を確認（正しく63になるはず）
length(unique_distributions)

# 結果を表示
print(unique_distributions)

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

763: １３２人目の素数さん [sage] 2025/03/25(火) 12:07:55.54 ID:JvJVqzfB

"
赤玉4個、黒玉3個、白玉2個、青玉1個を空箱なしで3つの箱にいれる。
箱を区別しないとき、入れ方は何通りあるか？
"
rm(list=ls())
# vector から remove　を除いた残りを返す
fn_rest <- function(vector,remove) Reduce(function(x, y) x[-match(y, x)], remove, init = vector)
fn_rest(c(1,2,2,3,3,3),c(2,3,3))

library(RcppAlgos)
v=rep(1:4,4:1)
pa=partitionsGeneral(10,3,repetition=TRUE) ; pa

pa[1,] # 1 1 8
c1=unique(comboGeneral(v,2)) ; c1
matrix(c("赤","黒","白","青")[c1],ncol=2) |> noquote()
re1=nrow(c1) ; re1

pa[2,] # 1 2 7
c2=unique(permuteGeneral(v,3)) ; c2
ans2=c2[apply(c2,1,\(x) x[2]<=x[3]),] ; ans2
matrix(c("赤","黒","白","青")[ans2],ncol=3) |> noquote()
re2=nrow(ans2) ; re2

pa[3,] # 1,3,6
c3=unique(permuteGeneral(v,4)) ; c3
ans3=unique(t(apply(c3,1,\(x) c(x[1],sort(x[2:4]))))) ; ans3
matrix(c("赤","黒","白","青")[ans3],ncol=4) |> noquote()
re3=nrow(ans3) ; re3

pa[4,] # 1 4,5
c4=unique(permuteGeneral(v,5)) ; c4
ans4=unique(t(apply(c4,1,\(x) c(x[1],sort(x[2:5]))))) ; ans4
matrix(c("赤","黒","白","青")[ans4],ncol=5) |> noquote()
re4=nrow(ans4) ; re4

pa[5,] # 2 2 6
c4=unique(comboGeneral(v,4)) ; c4
ans5=NULL
for(i in 1:nrow(c4)){
c22=c4[i,] ; c22
c2=unique(comboGeneral(c22,2)) ; c2
re22=cbind(c2,t(apply(c2,1, \(x) fn_rest(c22,x)))) ; re22
ans22=re22[apply(re22,1,\(x) sum(c(10,1)*x[1:2]) <= sum(c(10,1)*x[3:4])),] ; ans22
ans5=rbind(ans5,ans22)
}
matrix(c("赤","黒","白","青")[ans5],ncol=4) |> noquote()
re5=nrow(ans5) ; re5

pa[6,] # 2 3 5
c6=unique(permuteGeneral(v,5)) ; c6
ans6=unique(t(apply(c6,1,\(x) c(sort(x[1:2]),sort(x[3:5]))))) ; ans6
matrix(c("赤","黒","白","青")[ans6],ncol=5) |> noquote()
re6=nrow(ans6) ; re6

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

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

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

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

ぬこの手ぬこTOP 0.018s