高校数学の質問スレ(医者・東大卒専用) Part438 (991レス)
前次1-
抽出解除 レス栞

リロード規制です。10分ほどで解除するので、他のブラウザへ避難してください。
52: 132人目の素数さん [sage] 2024/08/15(木)19:23:48.55 ID:ysyyeRB1(2/2)
Rに移植

j=\(n){
a=sample(3,n,replace=TRUE)
count=1
while(length(unique(a))!=2){
a=sample(3,n,replace=TRUE)
count=count+1
}
b=sort(unique(a))
if(all(b==c(1,2))) winner=2
if(all(b==c(2,3))) winner=3
if(all(b==c(1,3))) winner=1
winners=sum(a==winner)
c(winners,count)
}

sim=\(n){
c=j(n)
winner=c[1]
count=c[2]
while(winner>1){
k=j(winner)
winner=k[1]
count=count+k[2]
}
count
}
res9=replicate(1e5,sim(9))
BEST::plotPost(res9,breaks='scott')
152: 132人目の素数さん [sage] 2024/08/25(日)12:59:59.55 ID:/qrXHaIo(3/11)
>>140
行列で計算させると算出時間が爆速なのにびっくりしました。
達人のスクリプトを改造して道具箱に保存しておきます。

(* j 人でジャンケンをしたときの終了までの回数の最頻値とその確率を返す *)
calc[j_]:=(
M=Table[If[n==m,1-(2^m-2)/3^(m-1),Binomial[m,n]/3^(m-1)],{m,1,j},{n,1,j}];
p=Differences@Table[MatrixPower[M,i][[j,1]],{i,0,10j}];
max=Max@p;
Flatten@{Position[p,max],max,N[max]}
)

In[3]:= calc[15]

Out[3]= {30, 65101358743766874914341259145354001254712997240185483777481087387163094008455949207206\

> 93198405902336999941273392239247993407703628004425130823658476122231689832344483758966236574\

> 375695 / 11947838420050013668726696739307151046843799152024135169583095938840977078626722578\

> 97327618239887790786549346048626664496721871548575328400043101228717425477619608889629973635\

> 327326175449, 0.0054488}
calc[15]
279: 132人目の素数さん [sage] 2024/10/31(木)13:47:25.55 ID:grkcalAP(3/9)
function (n, m)
{
stopifnot(is.numeric(n), is.numeric(m))
if (length(n) == 1) {
n <- rep(n, length(m))
}
else if (length(m) == 1) {
m <- rep(m, length(n))
}
ln <- length(n)
lm <- length(m)
if (ln != lm)
stop("Arguments 'n', 'm' must be scalars or have the same length.")
if (any(floor(n) != ceiling(n)) || any(floor(m) != ceiling(m)))
stop("Arguments 'n', 'm' must be integers or vectors of integers.")
k <- ifelse(m != 0, n%%m, NaN)
k <- ifelse(m != 0 & sign(n) != sign(m) & k != 0, k - m,
k)
return(k)
}
373
(1): 132人目の素数さん [sage] 2024/12/04(水)16:35:44.55 ID:M2J+bJMI(2/2)
Wolfram Language 14.0.0 Engine for Microsoft Windows (64-bit)
Copyright 1988-2023 Wolfram Research, Inc.

In[1]:= solve[nn_,n_,m_]:=(
ass=Counts[#]& /@ Flatten[Table[Sort /@ IntegerPartitions[x,{n},Range[0,m-1]],{x,m*Range[0,n-1]}],1];
li=KeyValueMap[List,#]& /@ ass;
{q,r}=QuotientRemainder[nn,m];
c=Flatten@{Table[q+1,r],Table[q,m-r]};
m2c[x_] :=(
i=If[x[[1]]==0,m,x[[1]]];
Binomial[c[[i]],x[[2]]]);
cases=Total[Times@@@ (m2c /@ # & /@ li)]; (* cases=Total[Times@@ #&/@ (m2c /@ # & /@ li)];*)
cases/Binomial[nn,n]
)

In[2]:= solve[100,10,5]

4555344594
Out[2]= -----------
22776722969
427
(1): 132人目の素数さん [sage] 2024/12/19(木)15:01:20.55 ID:K/JbHbND(1)
東大卒なら算出できるんじゃないの?
Wolfram,Pyth9n,C,Rのどれかは理工系卒なら使えて当然。
Fランは違うのか?
729
(3): 132人目の素数さん [sage] 02/20(木)18:42:43.55 ID:wRRdyhfp(2/3)
>>728
答は出せたの?Fランくん?
813: 132人目の素数さん [sage] 04/30(水)08:13:09.55 ID:wedVH8wl(10/10)
コメントが長すぎて読みにくくなった。
818: 132人目の素数さん [sage] 05/01(木)11:18:32.55 ID:L1qIlz9/(4/4)
ニッチな値の探索処理が終了しないコード

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)
828: 132人目の素数さん [sage] 05/05(月)19:09:02.55 ID:LDY1RQtT(3/4)
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.
847
(1): 132人目の素数さん [sage] 05/25(日)04:22:37.55 ID:P4nhnL8B(1/6)
# dbeta(L,a,b) == dbbeta(U,a,b)
# Solve[L^(a-1)(1-L)^(b-1)==U^(a-1)(1-U)^(b-1), b]

L=1/7
U=1/5
credMass = 0.95
f = function(a) 1 + ((a - 1) * log(U / L)) / log((1 - L) / (1 - U))
g = function(a) pbeta(U,a,f(a)) - pbeta(L,a,f(a)) - credMass
(re=uniroot(g,c(1,1e5)))
a=re$root
b=f(a)
c(a,b)
curve(g(x),1,1.5*a,bty="l") ; abline(h=0,lty=3)
a/(a+b) # mean
(a-1)/(a-1+b-1) # mode

library(LearnBayes)
beta.select(list(x=1/7,p=0.025),list(x=1/5,p=0.975))
923: 132人目の素数さん [sage] 06/22(日)13:20:05.55 ID:AY7cZjkg(4/10)
>>922
AIは解答を出しておりません
前次1-
スレ情報 赤レス抽出 画像レス抽出 歴の未読スレ AAサムネイル

ぬこの手 ぬこTOP 0.040s