臨床統計もおもしろいですよ、その3
library(rjags) s1=22;s2=10;n1=158;n2=78 c(s1/n1,s2/n2,s1/n1-s2/n2) shape1=shape2=0.5 data=list(s1=s1,s2=s2,n1=n1,n2=n2,shape1=shape1,shape2=shape2) modelString=' model{ # data s1 ~ dbin(p1,n1) s2 ~ dbin(p2,n2) # model p1 ~ dbeta(shape1,shape2) p2 ~ dbeta(shape1,shape2) delta = p1 - p2 # priors pri.delta= pri.p1 - pri.p2 pri.p1 ~ dbeta(shape1,shape2) pri.p2 ~ dbeta(shape1,shape2) } ' writeLines(modelString,'tmp.txt') n.iter=1e4 ; thin=1 jagsModel=jags.model('tmp.txt',data=data,n.chains=4,n.adapt=n.iter/5) update(jagsModel) codaSamples=coda.samples(jagsModel,c('delta','pri.delta'),n.iter,thin) plot(codaSamples) js=as.data.frame(as.matrix(codaSamples)) BEST::plotPost(js$delta,compVal = 0) library(polspline) fit=logspline(js$delta) pri.fit=logspline(js$pri.delta) xlim=c(-1,1) ; ylim=c(0,9) plot(fit,xlim=xlim, ylim=ylim,bty='l',lty=2) par(new=T) plot(pri.fit,xlim=xlim,ylim=ylim) dpost=dlogspline(0,fit) dprio=dlogspline(0,pri.fit) points(0,dpost) points(0,dprio,pch=19) dpost/dprio 転売ヤーの コロナ薬屋さん コンプ薬屋さん それが次の転売品ですか? 高音ぼったくりで通販サイトを荒らして 恥ずかしくないのですか? 転売品仕入れの時は夜中から並んで野糞をしてるって本当ですか? by 都内Sラン女子高生 http://itest.5ch.net/egg/test/read.cgi/hosp/1588253497/296 コロナ薬屋さん コンプ薬屋さんに共感されてますよ やっぱりお仲間だったんですね キンモー by 都内Sラン女子高生 今日も朝から人生が暇だけの コロナ薬屋 コンプ薬屋が 他スレでお医者さんごっこだぁW ちゃんとコテハン付けないとダメでちよぉぉぉ 今思うのは反日野党が政権をとってたら 日本はとっくに壊滅していたであろうこと だって反日なんだもの by 都内Sラン女子高生 それで 私立薬学部卒の Fラン統計ジジイ あらため コロナ薬屋さん あなたと違って大学は旧帝受験予定ですが それはそうと昨夜は何を転売するため並んでたんですか? 野糞しながらぼったくり転売は恥ずかしくないのですか? by 都内Sラン女子高生 全てのパチンコ屋を営業停止にして コロナ患者の収容施設にしよう 3人の女子大生に18種類の体位を教えてtヶ月後に覚えていた体位の数は以下の通りである. ?は欠測値。 記憶保持率をexp(-α*t)+βのモデルに当て嵌めて 4人めの冬子を含めて欠測値を推定せよ。 t=1 2 4 7 12 21 35 59 99 200 治子 18 18 16 13 9 6 4 4 4 ? 奈津子 17 13 9 6 4 4 4 4 4 ? 亜希子 14 10 6 4 4 4 4 4 4 ? 冬子 ? ? ? ? ? ? ? ? ? ? それで 私立薬学部卒の Fラン統計ジジイ あらため コロナ薬屋さん お仕事は転売ヤーで間違いないですね 頭が悪いのに転売のお仕事の為に統計の勉強をされたんですね by 都内Sラン女子高生 統計上は PCR件数と 死亡数は 正の相関関係がある と なるほど PCR信者はバカなんですね by 都内Sラン女子高生 今日も朝から人生が暇しかない Fラン統計ジジイ あらため コロナ薬屋 が 他スレでお医者さんごっこだぁW ちゃんとコテハン付けないとダメでちよぉぉぉ 今思うのは反日野党が政権をとってたら 日本はとっくに壊滅していたであろうこと だって反日なんだもの by 都内Sラン女子高生 ネットを見ると、大学の9月入学を支持する高3生の泣き言、みたいなニュースがありますね 私の周りの高3生にはそんなことを支持する者は誰もいないのですが どんなFラン高の生徒の話なんでしょうね 受験は生モノ 勉強時間が足りないのや試験日に体調が悪くて不合格になるのは、本人の実力のうちでしょうに by 都内Sラン女子高生 えぇぇえぇぇぇっ Fラン統計ジジイ あらため コロナ薬屋さんって 反日パヨクなんですかぁぁあ? サイテーですね 下水のゾウリムシよりもはるかにサイッテーですね ダンゴムシのごとくゴロゴロ転がりながら 海より深く反省しなさい by 都内Sラン女子高生 >>130 わりと、面倒。 library(rjags) t <- c(1, 2, 4, 7, 12, 21, 35, 59, 99, 200) nt <- length(t) slist <- 1:4 ns <- length(slist) k <- matrix(c(18, 18, 16, 13, 9, 6, 4, 4, 4, NA, 17, 13, 9, 6, 4, 4, 4, 4, 4, NA, 14, 10, 6, 4, 4, 4, 4, 4, 4, NA, NA, NA, NA, NA,NA,NA,NA,NA,NA, NA), nrow=ns, ncol=nt, byrow=T) colnames(k)=paste0('t',as.character(t)) rownames(k)=paste0('subject',1:4) n <- 18 dataList=list(t=t,k=k,nt=nt,ns=ns,n=n) cat("model { for(i in 1:ns){ for(j in 1:nt){ k[i,j] ~ dbin(theta[i,j],n) predk[i,j] ~ dbin(theta[i,j],n) } } for(i in 1:ns){ for(j in 1:nt){ theta[i,j] <- min(1 - 1e-16, exp(-alpha[i]*t[j]) + beta[i]) } # min(1, ....) => Node inconsistent with parents -- ERROR } for(i in 1:ns){ alpha[i] ~ dbeta(shape1a,shape2a) # dnorm(mu.a,lambda.a)T(0,1) beta[i] ~ dbeta(shape1b,shape2b) # dnorm(mu.b,lambda.b)T(0,1) } shape1a ~ dt(0,pow(2.5,-2),1)T(0,) # mu.a ~ dbeta(1,1) shape2a ~ dt(0,pow(2.5,-2),1)T(0,) # mu.b ~ dbeta(1,1) shape1b ~ dt(0,pow(2.5,-2),1)T(0,) # lambda.a ~ dgamma(0.001,0.001)T(0.001,) shape2b ~ dt(0,pow(2.5,-2),1)T(0,) # lambda.b ~ dgamma(0.001,0.001)T(0.001,) } ", file='tmp.txt') pk4=NULL for(i in 1:10){ pk4=cbind(pk4,eval(str2lang(paste0("js$\'predk[4," ,i, "]\'")))) } colnames(pk4)=c(paste0('t0',1:9),'t10') boxplot(pk4) head(pk4) df4=as.data.frame(pk4) head(df4) predk4=tidyr::pivot_longer(df4,col=1:10) head(predk4) colnames(predk4)=c('time','items') library(ggplot2) g <- ggplot(predk4,aes(time,items)) g + geom_violin() g + geom_boxplot() reshape2 が tidyr::pivot_longer(df4,col=1:10) になっていた。 使っていないと、ggplot2も忘れてしまうなぁ。 g <- ggplot(predk4,aes(time,items)) g + geom_boxplot() https://i.imgur.com/krliw0d.png バイオリンプロットというよりナメクジプロットだな。 g + geom_violin() https://i.imgur.com/fYZHA53.png histの$densityを使って面積比で反映させてみた https://i.imgur.com/bIigw6R.png 標準のboxplotでも十分な気がする。 https://i.imgur.com/L210Sf8.png えぇぇえぇぇぇっ Fラン統計ジジイ あらため コロナ薬屋さん および コンプ薬屋さんって 裏口容疑者なんですかぁぁあ? だと思いましたぁぁあぁぁ 早く出頭するんですよ あなたたちのレスの全てから 頭の悪さを感じますからねえぇぇぇ by 都内Sラン女子高生 交絡因子を考慮しての記憶減衰の階層ベイズモデル デバッグできたので、事前分布をいろいろ変えてみて再現性を確認。 内視鏡休診で纏まった勉強する時間が取れて( ・∀・)イイ!! ##### psychophysical function with confounder cat(' model{ for(i in 1:nsubjs){ for(j in 1:nstim[i]){ z[i,j] ~ dbern(phi[i]) pi[i,j] ~ dunif(0,1) r[i,j] ~ dbin(theta[i,j],n[i,j]) theta[i,j] <- ifelse(z[i,j]==1,ilogit(alpha[i] + beta[i]*(x[i,j]-xmean[i])),pi[i,j]) } } for(i in 1:nsubjs){ phi[i] <- phi(probit) alpha[i] ~ dnorm(mua,lambdaa) beta[i] ~ dnorm(mub,lambdab) } probit ~ dnorm(muphi,lambdaphi) muphi ~ dnorm(0,0.001) mua ~ dnorm(0,0.001) mub ~ dnorm(0,0.001) lambdaphi ~ dt(0,pow(2.5,-2),1)T(0,) # dunif(0.3) lambdaa ~ dt(0,pow(2.5,-2),1)T(0,) # dunif(0,1000) lambdab ~ dt(0,pow(2.5,-2),1)T(0,) # dunif(0,1000) } ',file='tmp.txt') 開業医ではないあなたが 先に開業医板を荒らすから 私がここに降臨して差し上げてるんですよ それすら判らないくらい知能が低いんですかねぇ 女子高生にブチ切れる異常者のジジイって どんなホラー映画の地獄絵図ですかW by 都内Sラン女子高生 >>143 half cauchy にしても一様分布にしても事後分布はほとんど変化はないな。 モデルが優れているからだろう。 さて今日も学校は終わり。 お家に帰ってお勉強しましょう。 あ〜あ、Fラン統計ジジイは やっぱりくるくるぱーだったんですねぇ by 都内Sラン女子高生 # 相関係数の差の検定 n1=150 r1=0.29 n2=120 r2=0.5 pnorm(-abs(atanh(r1)-atanh(r2))/sqrt(1/(n1-3)+1/(n2-3))) 離散量のプロットは重なってしまうと相関がわかりにくくなる。 jitterつけると重なりは防げるな。 頻度に応じて面積比でプロット https://i.imgur.com/QVrj8hW.png n=10 ; r=500 x=sample(n,r,rep=T,prob=rbeta(n,0.5,0.5)) y=sample(n,r,rep=T,prob=rbeta(n,0.5,0.5)) plot(rep(0,24),pch=1:24,cex=2,axes=F,ann=F,ylim=c(0,0.1)) ;axis(1,at=1:24) xy=paste0(x,y) cex=numeric() for(i in 1:length(x)){ cex[i]=sum(xy[i]==xy) } layout(matrix(1:4,2,b=T)) par(mar=c(1,1,1,1)) plot(x,y,pch=16) plot(x,jitter(y),pch=16) plot(jitter(x),jitter(y),pch=16) plot(x,y,pch=16,cex=sqrt(cex)) とりあえず、完成。多分、もっと便利なパッケージがあるんだろうな。 graphics.off() oldpar=par() par(mar=c(1,1,1,1),oma=c(3,3,3,3),bty='l',ann=F,pch=16) plot(rep(0,24),pch=1:24,cex=2,axes=F,ann=F,ylim=c(0,0.1)) ;axis(1,at=1:24) n=10 ; r=250 x=sample(n,r,rep=T,prob=abs(rnorm(n))) y=sample(n,r,rep=T,prob=rbeta(n,0.5,0.5)) plot(x,y) xy=paste0(x,y) cex=numeric() for(i in 1:length(x)){ cex[i]=sum(xy[i]==xy) } layout(matrix(1:4,2,b=T)) plot(jitter(x),jitter(y)) plot(x,y,cex=sqrt(cex)) plot(x,y,cex=3,col=rgb(0.1,0.1,0.1,0.25)) plot(x,y,cex=sqrt(cex),col=rgb(0.1,0.1,0.1,0.5)) layout(1) par(oldpar) 学校の先生に習いましたが、統計的にはPCR検査はむやみにやらないのが正解なんですね by 都内Sラン女子高生 【国民が誤解】安倍晋三「相談・受診の目安=PCR実施の基準ではない。国民への周知が足りなかった」 [スタス★] https://asahi.5ch.net/test/read.cgi/newsplus/1589284722/ 世界の頭脳は実態把握には大量検査が必要とわかっているなぁ。 東大の児玉名誉教授も同じ意見だな。 https://www.covid19-yamanaka.com/cont4/14.html 検診がメインだったから、内視鏡バイトは休診。 俺は6割の休業補償だが、銭ゲバ開業医はパート医を首にしているんだろうなぁ。 今日も朝から人生が暇しかない Fラン統計ジジイ あらため コロナ薬屋 が 発狂して連投だぁぁぁW ちゃんとコテハン付けないとダメでちよぉぉぉ 今思うのは反日野党が政権をとってたら 日本はとっくに壊滅していたであろうこと だって反日なんだもの by 都内Sラン女子高生 >>156 内視鏡バイトは休診なので、6割の休業補償はあるけど。 あんたのところは、感染拡大阻止に休診にして職員には休業補償してんの? まとまった時間がとれるので、https://bayesmodels.com/ の英文テキストで独学中。 行き詰まったら数学版で相談。 相関係数->FisherのZ変換->ベイズファクター r=js$r #The Fisher z-transformation: r.F <- atanh(r) fit.posterior <- logspline(r.F) plot(fit.posterior) posterior <- dlogspline(atanh(0), fit.posterior) r.F.prior <- atanh(runif(n.iter,-1,1)) fit.prior <- logspline(r.F.prior) prior <- dlogspline(atanh(0), fit.prior) BF10 <- prior/posterior BF10 まあ、ベイズファクターの計算結果はさほど変わらんな。 > BEST::plotPost(js$r,compVal = 0,showCurve = F,col='lightgreen') > lines(density(js$r),col='green') > abline(h=0.5,col=4) > fit.post=logspline(js$r) > points(0,dlogspline(0,fit.post),pch=19,col='white',cex=1.5) > points(0,dlogspline(0,fit.post),cex=1.5) > points(0,dunif(0,-1,1),pch=19,cex=1.5) > curve(dlogspline(x,fit.post),add=T) > fit.F=logspline(atanh(r)) > curve(dlogspline(x,fit.F), col=2, add=T) > d.prio=dunif(0,-1,1) > d.post=dlogspline(0,fit.post) > d.post/d.prio [1] 2.917773 > dlogspline(atanh(0),fit.F)/d.prio [1] 2.911407 あんまり相関係数の信頼区間など気にしていなかった。 必要ならmcmcかbootstrapすれば出せるし、 corr.testのソースを sink('cor.text.txt') stats:::cor.test.default sink() で覗いてみた。 if (n > 3) { if (!missing(conf.level) && (length(conf.level) != 1 || !is.finite(conf.level) || conf.level < 0 || conf.level > 1)) stop("'conf.level' must be a single number between 0 and 1") conf.int <- TRUE z <- atanh(r) sigma <- 1/sqrt(n - 3) cint <- switch(alternative, less = c(-Inf, z + sigma * qnorm(conf.level)), greater = c(z - sigma * qnorm(conf.level), Inf), two.sided = z + c(-1, 1) * sigma * qnorm((1 + conf.level)/2)) cint <- tanh(cint) FisherのZ変換して正規分布近似して信頼区間を出して、それをZ変換の逆関数で出しているな。 cat(' "名","本人","家族含む" "鳩山由紀夫首相",144269,144269 "菅直人副総理兼国家戦略担当相",905,2232 "原口一博総務相",914,1220 "千葉景子法相",3523,3523 "岡田克也外相",3273,8641 "藤井裕久財務相",14356,20214 "川端達夫文部科学相",4024,5583 "長妻昭厚生労働相",0,891 "赤松広隆農相",4864,5934 "直嶋正行経済産業相",3333,3333 "前原誠司国土交通相",741,1441 "小沢鋭仁環境相",2089,4014 "北沢俊美防衛相",309,609 "平野博文官房長官",1195,1875 "中井洽国家公安委員長",1296,1296 "亀井静香金融・郵政担当相",9427,18745 "福島瑞穂消費者・少子化担当相",12734,25000 "仙谷由人行政刷新担当相",1968,3987 ', file='shisan.csv') X=read.csv("shisan.csv") X library(boot) md = boot(X, function(df,idx) median(df[idx,2]), R=1e4) plot(md) mn = boot(X, function(df,idx) mean(df[idx,2]), R=1e4) plot(md) boxplot(list(median=md$t,mean=mn$t),horizontal = TRUE, las=1) >>162 パッケージ bootをつかうと library(boot) df=data.frame(x,k) bt=boot(df,function(df,idx) cor(df[idx,1],df[idx,2]),R=1e5) plot(bt) HDInterval::hdi(bt$t)[1:2] BEST::plotPost(as.vector(bt$t),xlab=bquote(cor),compVal = 0) > HDInterval::hdi(bt$t)[1:2] [1] -0.06736828 0.31237556 https://i.imgur.com/vCeo1Yn.png >>156 アベノミクスでの経済成長率は震災があった民主の頃より低いという現実を認められないんだなぁ。 日本を破壊したいなら、安倍一択。 安倍の長期政権で中国も韓国も笑いが止まらないはず、 日本が衰退する政策ばかりを選択しているから。 野党が反感をもたれて安倍信者が増えることこそ日本破壊が進んで パヨクの思う壺 結局、これが核心部分 Rt =I t (number of new infections generated at time t) / Σ[s=1,t] I t-s * Ws ( = total infectiousness of infected individuals at time t) Ws : an infectivity profile given by a probability distribution ws, dependent on time since infection of the case, s, but independent of calendar time, t. E[I t] = Rt Σ[s=1,t] I t-s * Ws Σ[s=1,t] I t-s * Wsの部分が畳み込み積分で Ws ∝ serial interval # https://www.nejm.org/doi/full/10.1056/NEJMoa2001316 #--- incubation period --- # from Li et al NEJM 2020 # lognormal mean = 5.2 ln.par1 = 1.434065 ln.par2 = 0.6612 x=rlnorm(1e6,ln.par1,ln.par2) BEST::plotPost(x) curve(dlnorm(x,ln.par1,ln.par2),add=T) quantile(x,c(0.025,0.5,0.975)) " c.onset = c.infected + c.incubation d.onset = d.infected + d.incubation c.onset - d.onset = onset.delay c.infected - d.infected + c.incubation - d.incubation = delay c.infected - d.infected = onset.delay + d.incubation - c.incubation " onset.delay = 2 d.incubation = rlnorm(1e6,ln.par1,ln.par2) c.incubation = rlnorm(1e6,ln.par1,ln.par2) infection.delay = onset.delay + d.incubation - c.incubation BEST::plotPost(infection.delay,compVal = 0) stancode= " data{ real onset_delay; real ln_par1; real ln_par2; } parameters{ real <lower=0> d_incubation; real <lower=0> c_incubation; } transformed parameters{ real infection_delay = onset_delay + d_incubation - c_incubation; } model{ d_incubation ~ lognormal(ln_par1,ln_par2); c_incubation ~ lognormal(ln_par1,ln_par2); } " model=stan_model(model_code = stancode) fn.stan <- function(delay){ dataList=list(onset_delay=delay,ln_par1=ln.par1,ln_par2=ln.par2) fit=sampling(model,data=dataList) ms=rstan::extract(fit) mean(ms$infection_delay < 0) } fn.stan(2) onset_delays=0:20 y=sapply(onset_delays,fn.stan) plot(onset_delays,y, ylab='Pr[ Infected Later ])',axes=F) ; axis(1) 今日も朝から人生が暇しかない Fラン統計ジジイ あらため コロナ薬屋 が 発狂して連投だぁぁぁW ちゃんとコテハン付けないとダメでちよぉぉぉ 今思うのは反日野党が政権をとってたら 日本はとっくに壊滅していたであろうこと だって反日なんだもの by 都内Sラン女子高生 runJags=run.jags('TEMPmodel.txt',monitor=c('p1','p2','diff'), data=dataList,n.chains=4,sample=10000,burnin=4000) coda::gelman.plot(runJags) codaSamples = as.mcmc.list(runJags) 今日も朝から人生が暇しかない Fラン統計ジジイ あらため コロナ薬屋 が 発狂して連投だぁぁぁW 相変わらず臨床経験や社会経験がゼロなのが 丸わかりのレスですねぇぇぇぇぇ ちゃんとコテハン付けないとダメでちよぉぉぉ 今思うのは反日野党が政権をとってたら 日本はとっくに壊滅していたであろうこと だって反日なんだもの by 都内Sラン女子高生 今日のくるくるぱーのIDだああぁぁあ id:3+Rla4zY くるくるぱーが反日野党に肩入れして 自滅するのは自業自得だけど 善良な日本人は巻き込まれたくないですねえ by 都内Sラン女子高生 今日も朝から人生が暇しかない Fラン統計ジジイ あらため コロナ薬屋 が 発狂して連投だぁぁぁW 相変わらず臨床経験や社会経験がゼロなのが 丸わかりのレスですねぇぇぇぇぇ 今思うのは反日野党が政権をとってたら 日本はとっくに壊滅していたであろうこと だって反日なんだもの by 都内Sラン女子高生 すいません、すごい基本的なことなのですが教えてください。 現在いる実験室で、様々細胞と遺伝子変異させた細胞に抗がん剤をかけて死亡率をみてます。 抗がん剤の量も比較してるのです 1: 細胞aとbに対して同量の抗がん剤を使用し死亡率を見る場合 2: 同種の細胞に違う濃度の抗がん剤をかけて死亡率を比較する場合。 両者とも何検体かとって平均の比較をする場合、対応するt検定ですか?、それともf値をみた上で対応しないt検定になるのですか? 今日も朝から人生が暇しかない Fラン統計ジジイ あらため コロナ薬屋 が 発狂して連投だぁぁぁW 身の程を弁えない謎の上から目線で政治談批評 老害アルアルですねぇぇぇぇぇ 今思うのは反日野党が政権をとってたら 日本はとっくに壊滅していたであろうこと だって反日なんだもの by 都内Sラン女子高生 日本衰退を願う勢力は安倍一択。 観光立国という亡国政策で衰退途上国。 観光でしか生きていけない国になったときに中国からの出国禁止すれば日本は中国に平伏すしかなくなる。 >>174 対数変換して検定した方がいい場合もある。 テレビを見てるんですが やっぱりマスコミはクズですね 問題のあった新聞社の社長は直ちに謝罪して 新聞社を解体するべきですね 反日野党は日本の足を引っ張ることしかしてませんね 存在意義がないとしか言えません だって反日なんだもの by 都内Sラン女子高生 >>178 対数変換した場合、検定方法は何使えばいいですか? ある人物Dが新型コロナ肺炎に罹患したとする。行動調査によって発症前にキャバクラに行っており 接客したキャバ嬢Cが人物D発症の2日後に発症していたことがわかった。 キャバ嬢Cは人物Dから移されたと主張して1億円の賠償を求めている。 潜伏期間には幅がありキャバ嬢Cから移された可能性もあると主張してその確率を計算して賠償金を値切りたい。 いくら値切れるか計算せよ。 #--- incubation period --- # from Li et al NEJM 2020 # lognormal mean = 5.2 ln_par1 = 1.434065 ln_par2 = 0.6612 Gt <- function(delay){ C=rlnorm(1e6,ln_par1,ln_par2) D=rlnorm(1e6,ln_par1,ln_par2) mean(C-D > delay) } Gt(2) library(polspline) delay=2 c=rlnorm(1e5,ln_par1,ln_par2) d=rlnorm(1e5,ln_par1,ln_par2) hist(c-d, freq=F,col='white',breaks=100,ylim=c(0,0.11), xlim=c(-30,30),ann=F,axes=F) ; axis(1) fit=logspline(c-d) curve(dlogspline(x,fit),-30,30,ann=F,bty='n',add=T) segments(delay,0,delay,dlogspline(delay,fit),pch=19,col=2) curve(dlogspline(x,fit),delay,30,add=T,type='h',col=2) 1-plogspline(delay,fit) fn <- function(delay) 1- plogspline(delay,fit) curve(fn(x),0,14, bty='n' ,xlab='Delay', ylab='Probability') >>181 t検定は二つのことを比べることしかできないですよね? >>183 多重検定やりたいならTukeyとか色々ある。 >>184 図の点線が対照のmockで実線が遺伝子変化させた細胞、 横軸が加えた薬剤の濃度です。 ルシフェラーゼの比を見ています 同じ濃度で、モックと遺伝子変化させた群を比較させたい時は統計は何使えばいいですか? https://i.imgur.com/VXfJ8TD.jpg 好きなの使えば。 厳しいbonferroni補正したpairwise.t.testとか 回帰係数の比較ならコクランアーミテッジだったかな。 Rだとprop.trend.testでできたはず! 皆さんありがとうございます すいませんjmpかエクセル統計で計算できるのだと嬉しいです。 Rだと kruskal.testで多重比較で有意差確認してから 平均比較ならpairwise.t.test 比率比較ならpairwise.prop.test 補正法は指定できる。 エクセルのマクロは売られている。 使ったこともないけど。 https://bellcurve.jp/ex/function/kruskal.html 今日も朝から人生が暇しかない Fラン統計ジジイ あらため コロナ薬屋 が 発狂して連投だぁぁぁW それにしても反日野党や反日マスコミは日本の足を引っ張ることしかしてませんね 存在意義がないとしか言えません だって反日なんだもの by 都内Sラン女子高生 " >勤務医が院長のおれの悪口を妻に話してたそうだ >腹が立ったからクビにする を登場人物としてみた。 (問題) 新型コロナに勤務医が罹患。 勤務医が発症した翌日に院長が発症、その2日後に妻が発症した。 (1)妻が感染源である(最初に感染していた)確率を求めよ。 (2)感染順が妻→勤務医→院長の順である確率を求めよ。 #--- incubation period --- # from Li et al NEJM 2020 # lognormal mean = 5.21 sd=3.86 # ln.par1 = 1.434065 # ln.par2 = 0.6612 rm(list=ls()) ln_par1 = 1.434065 ln_par2 = 0.6612 Aincub=rlnorm(1e6,ln_par1,ln_par2) # 勤務医 Bincub=rlnorm(1e6,ln_par1,ln_par2) # 院長 Cincub=rlnorm(1e6,ln_par1,ln_par2) # 妻 " Ainfected=Aonset-Aincub Binfected=Bonset-Bincub=Aonset+1-Bincub Cinfected=Conset-Cincub=Bonset+2-Cincub=Aonset+3-Cincub " Aonset=0 Ainfected=Aonset-Aincub Binfected=Aonset+1-Bincub Cinfected=Aonset+3-Cincub ABCinfected = as.data.frame(cbind(Ainfected,Binfected,Cinfected)) boxplot(ABCinfected) fn1 <- function(x) min(x)==x['Cinfected'] mean(apply(ABCinfected,1,fn1)) fn2 <- function(x){ x['Cinfected'] < x['Ainfected'] & x['Ainfected'] < x['Binfected'] } mean(apply(ABCinfected,1,fn2)) まず全裸になり ( : ) ( ゜∀゜)ノ彡 <( ) ノωヽ 自分の尻を両手でバンバン叩きながら白目をむき 从 Д゚ ) て ( ヾ) )ヾ て < < 人__人__人__人__人__人__人__人__人__人__人 Σ て Σ びっくりするほどユートピア! て人__人_ Σ びっくりするほどユートピア! て ⌒Y⌒Y⌒Y) て Y⌒Y⌒Y⌒Y⌒Y⌒Y⌒Y⌒Y⌒Y⌒Y⌒Y⌒ _______ |__ ヽ(゜∀゜)ノ |\_〃´ ̄ ̄ ヽ..ヘ( )ミ | |\,.-〜´ ̄ ̄ ω > (∀゜ )ノ \|∫\ _,. - 、_,. - 、 \ ( ヘ) \ \______ _\< \ || ̄ ̄ ̄ ̄ ̄ ̄ ̄ | \||_______ | これを10分程続けると妙な脱力感に襲われ、解脱気分に浸れる。 x1,x2,...,xnの順に発症 その間隔はt1,t2,..,tn-1としたときにxiが感染源であった確率を算出するプログラムを作れ。 安倍は新コロナは中国発と認めたから春節ウェルカムは外患誘致。 早くアメリカと足並み揃えて中国に損害賠償請求して国民に赦しを乞うべき。 観光立国という亡国政策で日本は衰退途上国。 次の世代の日本人は実習生として中国に出稼ぎがで立場が逆転。 日本人なら過労死自殺しても経営者殺害はしないから酷使される。 観光でしか生きていけない国になったときに中国からの訪日禁止すれば日本は中国に平伏すしかなくなる。 尖閣どころか沖縄を寄越せとか言われても差し出すことになりそうだな。 >>195 モデルは簡単なのでStanやJagsを使わずに完成 WhoInfectedFirst <- function( # 発症間隔から感染源の確率を計算 t=c(1,2), # 発症間隔 k=1e5 # 乱数発生数 ){ ln_par1 = 1.434065 # 潜伏期間対数正規分布パラメータ ln_par2 = 0.6612 n=length(t)+1 # 発症人数 # 潜伏期間 x_incub = matrix(rep(NA,n*k),ncol=n) for(i in 1:n) x_incub[,i] = rlnorm(k,ln_par1,ln_par2) # 感染日(一人目の発症日:0) x_infected = matrix(rep(NA,n*k),ncol=n) x_infected[,1]= -x_incub[,1] for(i in 2:n) x_infected[,i] = sum(t[1:(i-1)]) - x_incub[,i] # i番目の発症者が感染源の確率 fi <- function(i){ mean(apply(x_infected,1,function(x) which.min(x)==i)) } data.frame(p=round(sapply(1:n,fi),2)) } WhoInfectedFirst(c(1,2)) WhoInfectedFirst(c(1,0,1,0,0)) # 翌日2人発症 翌々日3人発症 やったね たえちゃん Fラン統計ジジイ あらため コロナ薬屋は 今日もクルクルパーだよ >>197 忘れないようにStanで書いてみた。 ################# MCMC by stan ############ stancode=' data { int<lower=0> n; vector[n-1] t; real ln_par1; real ln_par2; } parameters { real<lower=0> x_incub[n]; } model { target += lognormal_lpdf(x_incub|ln_par1,ln_par2); } generated quantities{ real x_infected[n]; x_infected[1]= -x_incub[1]; for(i in 2:n){ x_infected[i] = sum(t[1:(i-1)]) - x_incub[i]; } } ' stanmodel=stan_model(model_code = stancode) saveRDS(stanmodel,'cavaret.rds') t=c(1,0,1,0,0) n=length(t)+1 ln_par1 = 1.434065 ln_par2 = 0.6612 data=list(t=t,n=n,ln_par1=ln_par1,ln_par2=ln_par2) fit=sampling(stanmodel,data=data,iter=1e5,thin=1,chains=4, control=list(adapt_delta=0.95)) pars=c('x_infected') print(fit,pars=pars) stan_dens(fit,pars=pars,separate_chains = T) ms=rstan::extract(fit) head(ms$x_infected) x_infected=ms$x_infected fi <- function(i){ # i番目の発症者が感染源の確率 mean(apply(x_infected,1,function(x) which.min(x)==i)) } data.frame(p=round(sapply(1:n,fi),2)) > data.frame(p=round(sapply(1:n,fi),2)) p 1 0.26 2 0.18 3 0.18 4 0.13 5 0.13 6 0.13 JAGSだと library(rjags) cat(' model{ ln_par1 = 1.434065 # 潜伏期間対数正規分布パラメータ ln_par2 = 0.6612 # 潜伏期間 for(i in 1:n){ x_incub[i] ~ dlnorm(ln_par1,ln_par2) } # 感染日(一人目の発症日:0) x_infected[1]= -x_incub[1] for(i in 2:n){ x_infected[i] = sum(t[1:(i-1)]) - x_incub[i] } }',file='tmp.txt') t=c(1,0,1,0,0) n=length(t)+1 dataList=list(t=t,n=n) jagsModel=jags.model('tmp.txt',data=dataList, n.chains=4,n.adapt=5e3) update(jagsModel) codaSamples=coda.samples(jagsModel,n.iter=1e5,thin=1,var=c('x_infected')) x_infected=as.matrix(codaSamples) fi <- function(i){ # i番目の発症者が感染源の確率 mean(apply(x_infected,1,function(x) which.min(x)==i)) } data.frame(p=round(sapply(1:n,fi),2)) > data.frame(p=round(sapply(1:n,fi),2)) p 1 0.20 2 0.17 3 0.17 4 0.15 5 0.15 6 0.15 > WhoInfectedFirst(c(1,0,1,0,0)) # 翌日2人発症 翌々日3人発症 p 1 0.26 2 0.18 3 0.18 4 0.13 5 0.13 6 0.13 多数決でStanの値を採用しよう。 発症した順が必ずしも感染順ではないので、 職員が発症しても責めたりせずに悪いのは春節ウェルカムした安倍のせいと心をまとめよう。 # 武漢大量検査 fit=PCRs5a(6500000,218,SEN=0.6,SD1=0.2,SPC=0.999,SD2=0.001, iter=1e5,warmup=2e4,thin=10,chain=1)$fit print(fit,digits=7,prob=c(0.025,0.5,0.975),pars=c('spc','sen','prev','p')) stan_dens(fit, separate_chains = T) stan_plot(fit) stan_ac(fit) stan_trace(fit) ms=rstan::extract(fit) density2D(ms$sen,ms$spc) summary(ms$p*1e7) ; MODE(ms$p*1e7)[1] ; hdi(ms$p*1e7)[1:2] > print(fit,digits=7,prob=c(0.025,0.5,0.975),pars=c('spc','sen','prev','p')) Inference for Stan model: 451fbbe263cf71a87ea567fe0852d645. 1 chains, each with iter=100000; warmup=20000; thin=10; post-warmup draws per chain=8000, total post-warmup draws=8000. mean se_mean sd 2.5% 50% 97.5% n_eff Rhat spc 0.9999755 0.0000001 0.0000074 0.9999644 0.9999742 0.9999920 8082 1.00007 sen 0.4561154 0.0025058 0.2177783 0.0776906 0.4454970 0.8752057 7553 1.00002 prev 0.0000336 0.0000009 0.0000767 0.0000007 0.0000186 0.0001547 8059 1.00005 p 0.0000341 0.0000000 0.0000023 0.0000296 0.0000341 0.0000388 7840 0.99988 > summary(ms$p*1e7) ; MODE(ms$p*1e7)[1] ; hdi(ms$p*1e7)[1:2] Min. 1st Qu. Median Mean 3rd Qu. Max. 258 325 341 341 356 433 x 342.03 lower upper 296.90 388.22 検査は900万人らしいから最終的な陽性者数の95%信頼区間は267-349と出てきたな。 > summary(ms$p*9e6) ; MODE(ms$p*9e6)[1] ; hdi(ms$p*9e6)[1:2] Min. 1st Qu. Median Mean 3rd Qu. Max. 233 293 307 307 320 390 x 307.83 lower upper 267.21 349.40 >>182 乱数発生の代わりに数値積分で計算 gt = 2 # generation time (infection interval) mu = 1.434065 sg = 0.6612 # 差のpdfの公式 # pdf1 <- function(x) dlnorm(x,mu,sg) # pdf2 <- function(y) dlnorm(y,mu,sg) # ∫[-∞,∞] pdf1(x+y)*pdf2(y) dy f <- function(x,y) dlnorm(x+y,mu,sg)*dlnorm(y,mu,sg) vf=Vectorize(f,vectorize.args = 'y') pdf <- function(x) integrate(function(y) vf(x,y), rel.tol = 1e-14,-Inf,Inf)$value pdf=Vectorize(pdf) curve(pdf(x),-30,30) cdf <- function(x) integrate(pdf,x,Inf,rel.tol = 1e-14)$value cdf(gt) # 乱数発生との比較 c=rlnorm(1e6,mu,sg) d=rlnorm(1e6,mu,sg) hist(c-d, freq=F,breaks=200,ylim=c(0,0.11),col='skyblue', xlim=c(-30,30),ann=F,axes=F) ; axis(1) curve(pdf(x),add=T) mean(c-d > 2) gt = 2 # generation time (infection interval) mu = 1.434065 sg = 0.6612 # 差のpdfの公式 # pdf1 <- function(x) dlnorm(x,mu,sg) # pdf2 <- function(y) dlnorm(y,mu,sg) # ∫[-∞,∞] pdf1(x+y)*pdf2(y) dy f <- function(x,y) dlnorm(x+y,mu,sg)*dlnorm(y,mu,sg) vf=Vectorize(f,vectorize.args = 'y') pdf <- function(x) integrate(function(y) vf(x,y), rel.tol = 1e-14,-Inf,Inf)$value pdf=Vectorize(pdf) curve(pdf(x),-30,30) cdf <- function(x) integrate(pdf,x,Inf,rel.tol = 1e-14)$value cdf(gt) # 発症間隔を変化させてグラフ化 gts=seq(0,14,by=0.5) plot(gts,sapply(gts,cdf),type='l',xlab='発症間隔(日)', ylab='逆順感染確率',las=1) # c=rlnorm(1e6,mu,sg) d=rlnorm(1e6,mu,sg) hist(c-d, freq=F,breaks=200,ylim=c(0,0.11),col='skyblue', xlim=c(-30,30),ann=F,axes=F) ; axis(1) curve(pdf(x),add=T) mean(c-d > 2) >>210 ゴルゴ問題の答が出せるように勉強したのか? " 拳銃を一発撃ったときに、 狙った相手を撃ち殺す確率は、 Aは 1/3、Bは 1/2(50%)、Cは 1/1(100%)とします。 なお、この確率は、全員が知っているものとします。 拳銃を撃つ順番は、A、B、Cの順番で、 以降は最後の一人が生き残るまでこの順番を繰り返すものとします。 A、B、Cは拳銃を撃つときに誰を狙っても良いこととします。 ただし、一発で二人を狙うことはできません。 A、B、Cの生き残る確率を求めなさい。 " Duel <- function( pa=1/3, pb=1/2, pc=1/1, k=1e5 ){ qa=1-pa # probability of dead ; alive qb=1-pb qc=1-pc # simulation of 3 alive to 2 alive # each shoots at the superior sniper # Hence,A shoots C -> B shoots C -> C shoots B f3 <- function(){ a=b=c=1 bshoot=FALSE # is next shooter B? while(a+b+c==3){ c=rbinom(1,1,qa) # A shoots C if(c==0) bshoot=TRUE # if A kills C, next shooter is B c=rbinom(1,1,qb) # B shoots C if(c==1) b=rbinom(1,1,qc) # if C alive, C shoots B } return(list(abc=c(a,b,c),bshoot=bshoot)) } # simulation of 2 alive to 1 alive # Notice who has the right to shoot f2 <- function(){ re3=f3() # simulation to 2 suvivors a=re3$abc[1] # 1 1 b=re3$abc[2] # 1 0 c=re3$abc[3] # 0 1 while(a+b+c==2){ # while 2 alive if(c==0){ # when C dead if(re3$bshoot){ # when A killed C, B can shoot A a=rbinom(1,1,qb) if(a==0) return(c(a,b,c)) } else{ # when B killed C b=rbinom(1,1,qa) # A can shoot B if(b==0) return(c(a,b,c)) } } else{ # when C alive (C killed B, C shot at superior B) c=rbinom(1,1,qa) # A can shoot C if(c==0) return(c(a,b,c)) # if A killed C a=rbinom(1,1,qc) # when A missed C,C can shoot A if(a==0) return(c(a,b,c)) } } } re2=replicate(k,f2()) return(apply(re2,1,mean)) } Duel(1/3,1/2,1) Duel(0.3,0.4,0.5) # Debugged verision rm(list=ls()) # clear workspace p=c(1/3,1/2,1) # killing probablity q=1-p # survival probability # Three survivors to two survivors f32 <- function(init=1){ abc=c(1,1,1) # dead or alive n=length(abc) # initial survivors 3 s=sum(abc) # current survivors sh=init # initial shooter while(s==n){ # while 3 survivors tmp=p tmp[sh]=0 # set shooter for 0 probability (could be negative) target=which.max(tmp) # target : the superior sniper abc[target]=rbinom(1,1,q[sh]) # target dead(0) or alive(1) # next shooter : if next sniper alive next index, otherwise survived sniper sh=ifelse(abc[sh%%n+1]==1,sh%%n+1,(1:n)[-c(sh,target)]) s=sum(abc) # how many surviors left } list(abc=abc,sh=sh) } # demo f32() apply(replicate(1e5,f32()$abc),1,mean) # Two survivors to one survivor f21 <- function(){ # Three survivors to 2 survivor re32=f32() # two surviors & next shooter abc=re32$abc # two survivor sh=re32$sh # shooter sv=which(abc==1) # index of two suvivor target=sv[sv!=sh] # index of target s=sum(abc) while(s==2){ # while two survivors, mutual shooting abc[target]=rbinom(1,1,q[sh]) # target dead(0) or alive(1) s=sum(abc) # how many surviors left tmp=target # exhange shooter for target target=sh sh=tmp } abc } f21() apply(replicate(1e5,f21()),1,mean) rm(list=ls()) mu=0.58 # mean of reproductive number size=0.45 # dispersion parameter (prob = size/(size+mu)) # its probability Rt=rnbinom(1e5,size=size,mu=mu) # random numbers of negative binomial distribution hist(Rt,breaks = 'scott',freq=F,ann=F) # show its histgram sim <- function(n=10){ # simulation infectee=0 # initial value while(infectee!=n){ # while infectee is unequal to n infector=sample(n,1) # prior discrete uniform distribution of infector number infectee=sum(sample(Rt,infector)) # total number of infectee } return(infector) # when n infectee, return infector number } spreader=replicate(1e5,sim()) # simulation & calculation hist(spreader,freq=F,ylab='',axes=F,breaks='scott') ; axis(1) HDInterval::hdi(spreader)[1:2] # 95% credibility interval BEST::plotPost(spreader) # graph with 95%CI & mean quantile(spreader,c(0.025,0.5,0.975)) summary(spreader) sum(spreader==1)/length(spreader) # the probability of single super-spreader > quantile(spreader,c(0.025,0.5,0.975)) 2.5% 50% 97.5% 4 8 10 > summary(spreader) Min. 1st Qu. Median Mean 3rd Qu. Max. 1.000 7.000 8.000 8.039 10.000 10.000 > sum(spreader==1)/length(spreader) # the probability of single super-spreader [1] 0.0016 デバッグ版 R=0.58 # mean of reproductive number k=0.45 # dispersion parameter (prob = k/(k+R)) # its probability Rt=rnbinom(1e5,k,mu=R) # random numbers of negative binomial distribution hist(Rt,breaks = 'scott',freq=F,ann=F) # show its histgram sim <- function(n=10){ # simulation infected=0 # initial value while(infected!=n){ # while infected is unequal to n infector=sample(n,1) # prior discrete uniform distribution of infector number infectee=sum(sample(Rt,infector)) # number of infectee infected=infectee+infector # number of infected } return(infector) # when n infected, return infector number } spreader=replicate(1e5,sim()) # simulation & calculation hist(spreader,freq=F,ylab='',axes=F,breaks='scott') ; axis(1) HDInterval::hdi(spreader)[1:2] # 95% credibility interval BEST::plotPost(spreader) # graph with 95%CI & mean summary(spreader) sum(spreader==1)/length(spreader) # the probability of single super-spreader まあ、パーティー中に感染させる人数に再生産数を流用していいかは疑問ではあるが、 実効結果は https://i.imgur.com/ApDSW1J.png 95% CI > HDInterval::hdi(spreader)[1:2] # 95% credibility interval lower upper 3 9 中央値、平均値、四分位値 > summary(spreader) Min. 1st Qu. Median Mean 3rd Qu. Max. 1.000 6.000 7.000 6.844 8.000 10.000 1人のスーパースプレッダーの確率 > sum(spreader==1)/length(spreader) # the probability of single super-spreader [1] 0.00104 # n種類のガチャからm種類を集めるまでの期待値 collector <- function(n=100,m=80,print=TRUE){ library(gmp) x=(n-m+1):n x=as.bigq(x) y=sum(n/x) if(print) print(y) return(asNumeric(y)) } collector(5,4) collector(100,80) library(randtests) DEL <- permut(1:4,3,function(x) print(paste0(x,collapse = ''),quote=F)) まじで、エクセル、できればLibreOffice で出来る 医療統計の本を出して欲しい、この先生なら出来る。 リアル223からです、数式を記載されている方とは別人です。 >>223 エクセルのソルバーなんて使えたものじゃないね。 エクセルはグラフも貧弱だし。 医療統計の本なら朝倉から丹波某が沢山書いた本が沢山でているよ。 ベイズはWINBUGの話で古いけど。 " 容量8Lの袋と容量5Lの袋を使って池の水を丁度4L集めたい。 袋に目盛りはついていません。 袋から袋への移し替えは全量で行います。 池からとる水の量や池に捨てる水の量には制限はありません。 最初に片方に満たした作業を1回目として 丁度4Lを集めるのに最低何回の移動が必要か? " abura <- function( a7=8, b3=5, c5=4){ # starting from the bigger pitcher movea7 <- function(xy){ # start from c(a7,0) x=xy[1] ; y=xy[2] # x==a7 if(x==a7) re=c(a7-(b3-y),b3) # x==0 if(x==0) re=c(a7,y) # y==b3 if(y==b3) re=c(x,0) # y==0 if(y==0 & x!=a7){ if(x>=b3) re=c(x-b3,b3) else re=c(0,x) } return(re) } STATUS=status=c(a7,0) i=1 while(!identical(status,c(0,0))){# i=i+1 status=movea7(status) STATUS=rbind(STATUS,status) } rownames(STATUS)=1:nrow(STATUS) colnames(STATUS)=c(paste0(a7,'L'),paste0(b3,'L')) (Bigger=STATUS) # starting from the smaller pitcher moveb3 <- function(xy){ # start from c(0,b3) x=xy[1] ; y=xy[2] if(y==b3){ if(x<=(a7-b3)) re=c(x+b3,0) else re=c(a7, b3-(a7-x)) } if(y==0) re=c(x,b3) if(x==a7) re=c(0,y) if(x==0) re=c(y,0) return(re) } STATUS=status=c(0,b3) i=1 while(!identical(status,c(0,0))){ # stop at c(5,0) for solution i=i+1 status=moveb3(status) STATUS=rbind(STATUS,status) } rownames(STATUS)=1:nrow(STATUS) colnames(STATUS)=c(paste0(a7,'L'),paste0(b3,'L')) (Smaller=STATUS) min_Bigger=min(which(apply(Bigger,1,function(x) c5 %in% x))) min_Smaller=min(which(apply(Smaller,1,function(x) c5 %in% x))) list(Bigger=Bigger,Smaller=Smaller, min_Bigger=min_Bigger,min_Smaller=min_Smaller) } '%&%' <- function(x,y) paste0(x,y) > '√2 =' %&% sqrt(2) %&% ', π = ' %&% pi [1] "√2 =1.4142135623731, π = 3.14159265358979" random <- function(n,char=c(LETTERS,letters,0:9)){ re=paste0(sample(char,n),collapse='') cat(re,'\n') invisible(re) } > random() ujbDVK8pEtsNFTAo1BlkZP72cXRHrQvI0WdMmGyxOn3qz5awLhCg6S9eYJU4if onCone <- function(p,q){ # (p,q) any figure on cone x=sqrt(1-q^2)/pi*sin(pi*p/sqrt(1-q^2)) y=q z=sqrt(1-q^2)/pi - sqrt(1-q^2)/pi *cos(pi*p/sqrt(1-q^2)) + q list(x=x,y=y,z=z) } oncone <- function(p,q,α=pi/6){ # (p,q) 展開図上の座標、頂点の角度=2α PQ=sqrt(p^2+q^2) β=pi*sin(α) θ=atan(p/q) rdash=PQ*β/pi γ=PQ*θ/rdash x=rdash*sin(γ) y=PQ/(tan(α)*sqrt(1+tan(α)^-2)) z=rdash*(1-cos(γ)) c(x=x,y=y,z=z) } https://i.imgur.com/qOojo9N.png numlockキーをbackspaceキーに変更するレジストリ(要再起動) Windows Registry Editor Version 5.00 [HKEY_LOCAL_MACHINE\SYSTEM\CurrentControlSet\Control\Keyboard Layout] "Scancode Map"=hex:00,00,00,00,00,00,00,00,02,00,00,00,0e,00,45,00,00,00,00,00 ソースは https://www.souichi.club/technology/numlock-on/ oncone <- function(p,q,α=40*pi/180){ # (p,q) 展開図上の座標、頂点の角度=2α OA=sqrt(p^2+q^2) θ=Arg(p+1i*q) # == atan2(p+1i*q) β=pi*sin(α) # π*r=R*β ; r=R*sin(α) rdash=OA*sin(α) # r'== OA*β/π, r'*π == OA*β : hemi-circle==arc(B-c) δ=(pi/2-θ)/sin(α) # r'*δ == OA*(π/2-θ) ; δ == OA*(π/2-θ)/r' x=rdash*sin(δ) # Re(E) y=rdash/tan(α) # tan(α)=r'/y z=rdash-rdash*cos(δ) # IM(E-B) ; rdash+rdash(cos(π-δ)) c(x=x,y=y,z=z) } https://i.imgur.com/Nq1xRAc.png #原子数=放射能/崩壊定数=放射能/(log2/半減期)=放射能*半減期*log2 #放射能が同じなら原子数∝半減期 cesium_now <- function(Date=NULL,RCs134=1,RCs137=1,TCs134 = 2.0652,TCs137 = 30.16171){ t=ifelse(is.null(Date),as.numeric((Sys.Date()-as.Date("2011/3/11"))/365.2425) ,as.numeric(as.Date(Date)-as.Date("2011/3/11"))/365.2425) # mol ratio cs <- RCs134*TCs134*(1/2)^(t/TCs134) + RCs134*TCs137*(1/2)^(t/TCs137) cs0 <- RCs134*TCs134 + RCs134*TCs137 ratio=cs/cs0 # radioavtivity ratio # decay constant = log(2)/half-life # ratioactivity ∝ decay constant * mol CS <- (1/2)^(t/TCs134) + (1/2)^(t/TCs137) CS0 <- TCs134 + TCs137 Ratio=CS/CS0 list(year=t,mol_ratio=ratio,radioactivity_ratio=Ratio) } cesium_now() cesium_now(NULL,1,1,2,30) # 放射能比=1:1 TCs134 <- 2.0652 # 半減期(年) TCs137 <- 30.16171 CS <- function(t) (1/2)^(t/TCs134) + (1/2)^(t/TCs137) curve(CS(x),0,30) uniroot(function(t,u0=1/2) CS(t)/CS(0)- u0, c(0,30))$root # 問題「4/5より大きく5/6より小さい分数で、分母がいちばん小さい分数はなに?」 fn <- function(lo=4/5,up=5/6){ i=1 flg=FALSE while(flg==FALSE){ for(j in 1:ceiling(i*up)){ flg = lo<j/i & j/i<up if(flg==TRUE){ ans=paste0(j,'/',i) break } } i=i+1 } cat(ans) invisible(c(j,i)) } fn() fn(15/17,177/200) # πの近似分数 fn(355/113,22/7) 電話番号0531454900の詳細情報 - 電話番号検索 http://web.archive.org/web/20191110142958/https ://www.jpnumber.com/numberinfo_0531_45_4900.html https://i.imgur.com/A2YWEIJ.png VGB=integral(function(x) pi*green(x)^2, -r,b0) # volume green&blue VR=integral(function(x) pi*red(x)^2, b0,a+R) # volume red AG=integral(function(x) 2*pi*green(x), -r,b0) # area green AB=integral(function(x) 2*pi*blue(x),b0,c+s) # area blue AR=integral(function(x) 2*pi*red(x)^2, c+s,a+R) # area red >>239 では、同業者の意見を聞いてみましょう。 https://egg.5ch.net/test/read.cgi/hosp/1592662437/73 73 卵の名無しさん sage 2020/06/23(火) 13:24:47.79 ID:riQXI/fH 宮廷卒だけど、一括りに医師免許と言ってるが、私大卒など医者とは思ってへんよ 私大入学というインチキを経由したイシャモドキが、あんま調子のんなや 面白いねたしかに AIみたいに定型文かコピペしか送ってこないんだもん じゃあそんな君にこれをプレゼント! うりゅうひろゆきの本スレにようこそ! このウンコスレは基本的にさらしアゲて使用します! ココは医師真似事務員ジイさんが、医師妬みの挙句、テメェで勝手に自称で医科歯科大を卒業した医師って設定で 妄想願望日記を発表するスレです お医者さんのマトモなレスなど全く不要! ゴミ箱まがいの公衆便所掃き溜め痰壷下品スレとして使用しましょうね そんなわけで、いまだ懲りない恥さらし(自称)医科歯科事務員ジジイ登場〜〜 │ 彡川川川三三三ミ〜 プウゥ〜ン | 川|川/ \|〜 ポワ〜ン ________ | ‖|‖ ◎---◎|〜 / | 川川‖ 3 ヽ〜 < 僕はうりゅうひろゆき | 川川 ∴)д(∴)〜 \________ | 川川 〜 /〜 カタカタカタ | 川川‖ 〜 /‖ _____ | 川川川川___/‖ | | ̄ ̄\ \ | / \__| | | ̄ ̄| | / \医科歯科 | | |__| | | \ |つ |__|__/ / | / 医師うらやましいぜ・・ | ̄ ̄ ̄ ̄| 〔 ̄ ̄〕 よし、医科歯科の医者と称して 私立のお医者様を攻撃だ!なんたって、医科歯科の看板借りてんだからな え?証明???自称だから医科歯科の学生証、卒業証書あるわけねえよ 同期の医師あげてみろ、ったて、そもそも医師じゃねえし居ねえよバカ でも病院でもらった、医科歯科の封筒はもってま〜す!これを証拠にしてやるぜ 専門医?欲しいけど・・いいや!そんなの必要ねえ! 俺は事務員だし 某院でn人の患者感染。死亡率pとすると 予想される死亡者数とその95%信頼区間はいくらになるか? " p=0.20 n=100 (CI=HDInterval::hdi(qbinom,size=n,prob=p)) y=dbinom(0:n,n,p) plot(0:n,y,type='h',bty='l') sum(dbinom(CI[1]:CI[2],n,p)) ウリュウのジジイ、5chでも社会でもまるで相手にされず! 1 meter * 1 meter * 1 milli meter (water) = 1kg 1 mmHg柱/m2 = 13.6kg重/m2 13.6*9.8 N/m2 = 133.3 Pa =1.33hPa 1 hPa = 1/1.33 mmHg柱 = 0.752 mmHg柱 1 mmH20 = 1 kg重/m2 1*9.8 N/m2 = 9.8 Pa = 0.98 hPa 1 hPa = 1/0.98 = 1.02 mmH2O https://news.mynavi.jp/article/20160317-a313/ about:config privatebrowsing.autostartと入力して設定対象を絞り込み ほんと庭なしバカ安達はくそ犬はうるせーわ道路でばかでかい声で会話するわやりたい放題 家族みんなデブだし 【HKEY_CURRENT_USER > Control Panel > Personazation > Desktop Slideshow】 tryCatch has a slightly complex syntax structure. However, once we understand the 4 parts which constitute a complete tryCatch call as shown below, it becomes easy to remember: expr: [Required] R code(s) to be evaluated error : [Optional] What should run if an error occured while evaluating the codes in expr warning : [Optional] What should run if a warning occured while evaluating the codes in expr finally : [Optional] What should run just before quitting the tryCatch call, irrespective of if expr ran successfully, with an error, or with a warning tryCatch( expr = { # Your code... # goes here... # ... }, error = function(e){ # (Optional) # Do this if an error is caught... }, warning = function(w){ # (Optional) # Do this if an warning is caught... }, finally = { # (Optional) # Do this at the end before quitting the tryCatch structure... } ) log_calculator <- function(x){ tryCatch( expr = { message(log(x)) message("Successfully executed the log(x) call.") }, error = function(e){ message('Caught an error!') print(e) }, warning = function(w){ message('Caught an warning!') print(w) }, finally = { message('All done, quitting.') } ) } Windows Registry Editor Version 5.00 [HKEY_CURRENT_USER\Control Panel\Desktop\WindowMetrics] 古い低スペックPCに、Insider Programの正規版でWindows 11を入れてみました。 ↓ Windows Insider Program公式サイト https://insider.windows.com/ja-jp/ <レジストリの変更> HKEY_LOCAL_MACHINE\SOFTWARE\Microsoft\WindowsSelfHost\UI\Selection UIBranch ← Dev UIContentType ←Mainline UIRing ← External HKEY_LOCAL_MACHINE\SOFTWARE\Microsoft\WindowsSelfHost\Applicability BranchName ← Dev ContentType ←Mainline Ring ← External HKEY_LOCAL_MACHINE\SYSTEM\Setup 新しいキー LabConfig を追加 LabConfigに以下のDWORD(32ビット値)を追加 BypassTPMCheck ← 1 BypassSecureBootCheck ← 1 Windows Updateでエラー0cx190012eが発生した場合、 C:\Users\Default\AppData\Local\Microsoft\Windows\WSUS からsetupconfig.iniを削除します。 AppraiserRes.dllのダウンロードサイト https://github.com/CodeProf14/Fix-TPM C:\$WINDOWS.~BT\SourcesのDLLと置き換える '%IN%' <- function(v,M) any(rowSums(M == v[col(M)]) == ncol(M)) c(4,5,6) %IN% matrix(c(1,2,3, 4,5,6, 7,8,9), nrow=3,ncol=3,byrow=TRUE) >>1 【全スレに注意連絡】 医療空間で連携して悪質な行為をする創.価.学.会.員.に注意してください。 【創.価.学.会.員.が連携して特定の患者をいじめたり、嫌がらせ、ガスライティング、不審な診察や治療、病院に受診できなくさせたり、たらい回し】 そのようなことが絶対にないように、みんなで不審な学.会.員の医療従事者の行動を見守りましょう。 そしておかしなことがあれば、そのまま放置せず、各所に報告や通報、ネットなどにも声をあげていきましょう。 隠蔽されやすい連携した嫌がらせは、そのように表沙汰にすることで弱体化していきます。 みんなで声をあげて悪質な連携を阻止していきましょう。 関連スレッド↓ 医療法人アイクレセント いなだ小児科・アレルギー科 【佐賀】2 https://itest.5ch.net/egg/test/read.cgi/hosp/1624113085 集団ストーカー自称おおせき https://lavender.5ch.net/test/read.cgi/minor/1622387124/l50 佐賀大学医学部附属病院 https://egg.5ch.net/test/read.cgi/hosp/1622027754/l50 情報交換などはこのスレでもできます↓ ◆◆◆◆◆◆◆ 医療空間での創.価.学.会 ◆◆◆◆◆◆◆ https://egg.5ch.net/test/read.cgi/hosp/1628326388/l50 情報のリークはこのスレにも出来ます↓ 誰にも言えない事を吐くスレ https://egg.5ch.net/test/read.cgi/hosp/1627796592/l50 hdkdl [R] retrieving object names passed indirectly to a function https://stat.ethz.ch/pipermail/r-help/2010-September/252840.html foo <- function(p1, p2) { p1name <- deparse(substitute(p1)) p2name <- deparse(substitute(p2)) # The real function runs a simulation, but # here we just confirm the object names # and values passed as args cat("passed", p1name, "with value", p1, "\n") cat("passed", p2name, "with value", p2, "\n") } " レミフェンタニルの添付文書には,BMIが25を超える肥満患者には,標準体重を用いた投与量の調整を行うことが推奨されている. 高齢者:血圧低下などの副作用があらわれやすいため,開始用量の半減など, 患者の全身状態を観察しながら投与量を決定 " # 2mg/20mL lwr=0.25μg/kg/min upr=0.05/kg/min ideal BMI=22 Ultiva<-function(cm,kg,age,lwr=0.25,upr=0.50,conc=2000/20,bmi=22){ BMI=kg/(cm/100)^2 LBM=ifelse(BMI>25,bmi*(cm/100)^2,kg) # lean body mass L=lwr*LBM*60/conc U=upr*LBM*60/conc Lower=ifelse(age<75,L,L/2) # mL/h Upper=ifelse(age<75,U,U/2) round(c(BMI=BMI,LBM=LBM,'L(mL/h)'=Lower,'U(mL/h)'=Upper),2) } Ultiva(167.5,77.5,33) Ultiva(137.2,50.8,87) " 持続投与量の目安は、セボフルランでは3-4 μg/kg/min、プロポフォールでは7.5 μg/kg/minと言われる。 添付文書:持続注入により投与する場合は、7μg/kg/分の投与速度で持続注入を開始する。 " Eslax<-function(kg,sevo=FALSE,lwr=0.6,upr=0.9,conc=50e3/5){ # 50mg/5mL precurarization=kg*0.03*1e3/conc l=lwr*1e3*kg/conc # mL bolus u=upr*1e3*kg/conc L=ifelse(sevo,3*kg*60/conc,7*kg*60/conc) # on label U=ifelse(sevo,4*kg*60/conc,7.5*kg*60/conc) # off label cat('precararization(mL) =',precurarization,'\n') cat('bolus(mL) =',l,'-',u,'\n') cat('continuous(mL/h) =',L,'-',U,'\n') } Eslax(77.5,sevo=FALSE) Eslax(50.8,sevo=TRUE) # Practically loading kg(mL/h) for 10 min and maintenance kg/10(ml/h) Precedex<-function(kg,conc=200/50,lwr=0.2,upr=0.7,load=6){ I=kg*load/conc # 初期(10分)負荷(添付文書) L=kg*lwr/conc # 維持下限(添付文書) U=kg*upr/conc # 維持上限(添付文書) i=kg*4/conc # 推奨初期負荷 m=kg*0.4/conc # 推奨維持 d=kg*0.1/conc # 譫妄予防 cat('sedation(mL/h) =',round(L,2),'-',round(U,2),'\n') cat('delirium prevention(mL/h) =',round(d,2),'\n\n') invisible(list(初期負荷=I,維持下限=L,維持上限=U,推奨初期負荷=i,推奨維持=m,譫妄予防=d)) } Precedex(50) > Precedex(50) sedation(mL/h) = 2.5 - 8.75 delirium prevention(mL/h) = 1.25 # pdfからcdfの逆関数を作ってhdiを表示させて逆関数を返す # 両端での演算を回避 ∫[0,1]は∫[1/nxxx,1-1/nxx]で計算 pdf2hdi <- function(pdf,xMIN=0,xMAX=1,cred=0.95,Print=TRUE,nxx=1001){ xx=seq(xMIN,xMAX,length=nxx) xx=xx[-nxx] xx=xx[-1] xmin=xx[1] xmax=xx[nxx-2] AUC=integrate(pdf,xmin,xmax)$value PDF=function(x)pdf(x)/AUC cdf <- function(x) integrate(PDF,xmin,x)$value ICDF <- function(x) uniroot(function(y) cdf(y)-x,c(xmin,xmax))$root hdi=HDInterval::hdi(ICDF,credMass=cred) print(c(hdi[1],hdi[2]),digits=5) if(Print){ par(mfrow=c(3,1)) plot(xx,sapply(xx,PDF),main='pdf',type='h',xlab='x',ylab='Density',col='lightgreen',bty='l') legend('top',bty='n',legend=paste('HDI:',round(hdi,3))) plot(xx,sapply(xx,cdf),main='cdf',type='h',xlab='x',ylab='Probability',col='lightblue',bty='l') pp=seq(0,1,length=nxx) pp=pp[-nxx] pp=pp[-1] plot(pp,sapply(pp,ICDF),type='l',xlab='p',ylab='x',main='ICDF',bty='l') par(mfrow=c(1,1)) } invisible(ICDF) } # ベータ分布の差の分布 beta.diff<-function(r1,r2,n1,n2,a=1,b=1){ f <- function(x,y) dbeta(x+y,a+r1,b+n1-r1)*dbeta(y,a+r2,b+n2-r2) vf=Vectorize(f,vectorize.args = 'y') pdf <- function(x) integrate(function(y) vf(x,y),-1,1)$value pdf=Vectorize(pdf) invisible(pdf) } pdf=beta.diff(r1,r2,n1,n2) MAC2Sevo=function(MAC,age,N2O=0,MAC40sevo=1.8){ GOS=function( age=85, # 年齢 FE_sevo=0.5, # 呼気セボフルラン濃度(%) FE_N20=66.6){ # 呼気笑気濃度(%) MAC40_sevo=MAC40sevo # 40歳セボフルランMAC(1.71,1.8,2.1と諸説あり) MAC40_N2O=104 # 笑気MAC MACage_sevo=MAC40_sevo*10^(-0.00269*(age-40)) MACage_N2O=MAC40_N2O*10^(-0.00269*(age-40)) MACage_total = FE_sevo/MACage_sevo + FE_N20/MACage_N2O return(MACage_total) # 混合ガスのMAC } f=function(x) GOS(age,x,N2O) - MAC f=Vectorize(f) c('FEsevo(%)'=round(uniroot(f,c(0,5))$root,2)) } MAC2Sevo(1.5,57) # 年齢を設定してMACを笑気とセボフルレン濃度でMACをグラフ化 Age=87 GOS=function( age=85, # 年齢 FE_sevo=0.5, # 呼気セボフルラン濃度(%) FE_N20=66.6, # 呼気笑気濃度(%) MAC40sevo=1.8){ # 40歳セボフルランMAC(1.71,1.8,2.1と諸説あり) MAC40_sevo=MAC40sevo MAC40_N2O=104 # 笑気MAC MACage_sevo=MAC40_sevo*10^(-0.00269*(age-40)) MACage_N2O=MAC40_N2O*10^(-0.00269*(age-40)) MACage_total = FE_sevo/MACage_sevo + FE_N20/MACage_N2O return(MACage_total) # 混合ガスのMAC } GOS=Vectorize(GOS) sev=seq(0,3,by=0.1) N2O=seq(0,70,by=1) MAC=outer(sev,N2O,function(x,y)GOS(Age,x,y)) contour(sev,N2O,MAC,nlevels=50, bty='l',xlab='Sevoflurane',ylab='nitrous oxide',main='MAC for age 87') contour(sev,N2O,MAC,level=0.75,lwd=2,add=T) f=\(n){ dat=c(0,rep(100,n-1)) m=mean(dat) s=sd(dat) 50+(0-m)/s } nn=1000:5000 plot(nn,sapply(nn,f)) uniroot(f,c(1000,10000)) f(2501) f(2502) calc = \(L,m,n){ f1 = \(x){ y = as.character(x) strsplit(y,'') |> unlist() |> as.numeric() } f2 = \(x){ y = f1(x) x%%n==0 & all(y==sort(y)) & sum(y)==m } x=((10^L-1)/9):(10^L-1) x[sapply(x,f2)] } calc(3,16,2) calc(5,24,3) FF 83 F8 03 72 FF 83 F8 03 EB pH=7.32 pCO2=23 HCO3=16 pH_l=7.36 pH_u=7.44 pCO2_l=36 pCO2_u=44 HCO3_l=21 HCO3_u=27 if(pH<pH_l & pCO2<pCO2_l){ primary="metabolic acidosis" expected_pCO2=40 + 1.2*(HCO3-24) if(pCO2<expected_pCO2) additional= "secondary respiratory alkalosis") if(pCO2>expected_pCO2) additional= "secondary respiratory acidosis") } if(pH>pH_u & pCO2>pCO2_u){ primary="metabolic alkalosis" expected_pCO2=40 + 0.7*(HCO3-24) if(pCO2<expected_pCO2) additional= "secondary respiratory alkalosis") if(pCO2>expected_pCO2) additional= "secondary respiratory acidosis") } if(pH<pH_l & pCO2>pCO2_u){ primary="respiratory aciosis" status=NULL if(HCO3_l<HCO3 & HCO3<HCO3_u) status="acute" else status="chronic" if(status=="acute") expected_HCO3=24+0.1*(pCO2-40) if(status=="chronic") expected_HCO3=24+0.4*(pCO2-40) if(HCO3<expectedHCO3) additional="secondary metabolic acidosis") if(HCO3>expectedHCO3) additional="secondary metabolic alkalosis") } if(pH>pH_u & pCO2<pCO2_l){ primary="respiratory alkalosis" status=NULL if(HCO3_l<HCO3 & HCO3<HCO3_u) status="acute" else status="chronic" if(status=="acute") expected_HCO3=24+0.2*(40-pCO2) if(status=="chronic") expected_HCO3=24+0.4*(40-pCO2) if(HCO3<expectedHCO3) additional="secondary metabolic acidosis") if(HCO3>expectedHCO3) additional="secondary metabolic alkalosis") cat(paste(status,primary),'\n') cat(additional) if(pH<pH_l & pCO2>pCO2_u){ primary="respiratory aciosis" status=NULL if(HCO3_l<HCO3 & HCO3<HCO3_u) status="acute" else status="chronic" if(status=="acute") expected_HCO3=24+0.1*(pCO2-40) if(status=="chronic") expected_HCO3=24+0.4*(pCO2-40) if(HCO3<expectedHCO3) additional="secondary metabolic acidosis") if(HCO3>expectedHCO3) additional="secondary metabolic alkalosis") } if(pH>pH_u & pCO2<pCO2_l){ primary="respiratory alkalosis" status=NULL if(HCO3_l<HCO3 & HCO3<HCO3_u) status="acute" else status="chronic" if(status=="acute") expected_HCO3=24+0.2*(40-pCO2) if(status=="chronic") expected_HCO3=24+0.4*(40-pCO2) if(HCO3<expectedHCO3) additional="secondary metabolic acidosis") if(HCO3>expectedHCO3) additional="secondary metabolic alkalosis") cat(paste(status,primary),'\n') cat(additional) Astrap <- \(pH,pCO2,HCO3){ pH_l=7.36 pH_u=7.44 pCO2_l=36 pCO2_u=44 HCO3_l=22 HCO3_u=26 if(pH_l<pH & pH<pH_u){ if(pCO2>pCO2_u){ primary="mixed respiratory acidosis" additional="and metabolic alkalosis" } if(pCO2<pCO2_l){ primary="mixed respiratory alkalosis" additional="and metabolic acidosis" } } if(pCO2_l<pCO2 & pCO2<pCO2_u){ if(pH>pH_u){ primary="mixed metaboic alkalosis" additional="and respiratory acidosis" } if(pH<pH_l){ primary="mixed metabolic acidosis" additional="and respiratory alkalosis" } } if(pH<pH_l & pCO2<pCO2_l){ primary="primary metabolic acidosis" expected_pCO2=40 + 1.2*(HCO3-24) if(pCO2<expected_pCO2) additional= "secondary respiratory alkalosis" if(pCO2>expected_pCO2) additional= "secondary respiratory acidosis" } if(pH>pH_u & pCO2>pCO2_u){ primary="primary metabolic alkalosis" expected_pCO2=40 + 0.7*(HCO3-24) if(pCO2<expected_pCO2) additional= "secondary respiratory alkalosis" if(pCO2>expected_pCO2) additional= "secondary respiratory acidosis" } status=NULL if(pH<pH_l & pCO2>pCO2_u){ primary="primary respiratory aciosis" if(HCO3_l<HCO3 & HCO3<HCO3_u) status="acute" else status="chronic" if(status=="acute") expected_HCO3=24+0.1*(pCO2-40) if(status=="chronic") expected_HCO3=24+0.4*(pCO2-40) if(HCO3<expected_HCO3) additional="secondary metabolic acidosis" if(HCO3>expected_HCO3) additional="secondary metabolic alkalosis" } if(pH>pH_u & pCO2<pCO2_l){ primary="primary respiratory alkalosis" if(HCO3_l<HCO3 & HCO3<HCO3_u) status="acute" else status="chronic" if(status=="acute") expected_HCO3=24+0.2*(40-pCO2) if(status=="chronic") expected_HCO3=24+0.4*(40-pCO2) if(HCO3<expected_HCO3) additional="secondary metabolic acidosis" if(HCO3>expected_HCO3) additional="secondary metabolic alkalosis" } cat(status,primary,'\n') cat(additional) } > Astrap(pH=7.23,pCO2=23,HCO3=16) primary metabolic acidosis secondary respiratory alkalosis > Astrap(pH=7.54,pCO2=23,HCO3=17) chronic primary respiratory alkalosis secondary metabolic acidosis > Astrap(pH=7.40, pCO2=23, HCO3=20) mixed respiratory alkalosis and metabolic acidosis > Astrap(pH=7.40, pCO2=23, HCO3=35) mixed respiratory alkalosis and metabolic acidosis > Astrap(pH=7.40, pCO2=23, HCO3=17) mixed respiratory alkalosis and metabolic acidosis >>1 【全スレ注意連絡】※告発内容はコピー自由カクサン希望 【裏で集団ストーカー犯罪を指示してる医師、稲田成安に気をつけて】 http://www.inada-kids.com/wordpress/wp-content/uploads/greeting.jpg 佐賀県のいなだ小児科・アレルギー科 の小児科医、稲田成安(いなだ しげやす)に気をつけてください 自分に批判的な医療従事者や患者、その他の人たちを、探偵やソウカガッカイの集団ストーカーを使って調べ上げて、嫌がらせをやらせて仕返し、弾圧してきます あらゆる嫌がらせで自殺や自滅に追い込む工作を仕掛けてきます (※ 5ちゃんねるでは「ソウカガッカイ」という言葉を漢字四字で投稿するとアクセス禁止になる事があるのでカタカナ表記にしてあります) もし稲田による集団ストーカー犯罪にあったら、証拠を取り、周囲や各機関に声を上げていきましょう。 そうする事で集団ストーカー組織は弱体化していきます 詳細は、↓のスレ 【ソウカ在日】稲田成安(いなだ しげやす) https://egg.5ch.net/test/read.cgi/hosp/1644311065/ いなだ小児科・アレルギー科 【佐賀】 https://egg.5ch.net/test/read.cgi/hosp/1622365146/-100 黒猫世(くろねこよ) https://rio2016.5ch.net/test/read.cgi/police/1648800949/l50 gyhjoooo 「 いなだ小児科・アレルギー科 集団ストーカー犯罪 」で検索 rm(list=ls()) library(RcppAlgos) pm=permuteGeneral(v=1:3,freqs=c(3,3,4)) n=nrow(pm) f=\(x){ t1=(1:10)[x==1] t2=(1:10)[x==2] t3=(1:10)[x==3] if(t1[1]>t2[1]){ temp=t1 t1=t2 t2=temp } list(t1,t2,t3) c(paste(LETTERS[t1],collapse=''),paste(LETTERS[t2],collapse=''), paste(LETTERS[t3],collapse='')) } ans=unique(t(apply(pm,1,f))) print(head(ans,15),quote=FALSE) print(tail(ans,15),quote=FALSE) " 10人の子供(A,B,C,....,J)を3チームに割り当てる。 各チームは1人以上の子供が割り当てられる。 何通りの分け方があるか? 例. (A) (J) (B,C,D,E,F,G,H,I) (A,C) (H,I,J) (B,D,E,F,G) " rm(list=ls()) gr=expand.grid(1:8,1:8) gr1=gr[rowSums(gr)<10,] f=\(x) 10-sum(x) gr2=cbind(gr1,apply(gr1,1,f)) gr3=unique(t(apply(gr2,1,sort))) gr3 ans=apply(gr3,1,\(x) choose(10,x[1])*choose(10-x[1],x[2])/ length(unique(x)) ) sum(ans) ans "10人の子供(A,B,C,....,J)を3人、3人、4人のチームに割り当てる。 何通りの分け方があるか、列挙せよ" rm(list=ls()) library(RcppAlgos) pm=permuteGeneral(v=1:3,freqs=c(3,3,4)) f=\(x){ t1=(1:10)[x==1] t2=(1:10)[x==2] t3=(1:10)[x==3] if(t1[1]>t2[1]){ temp=t1 t1=t2 t2=temp } list(t1,t2,t3) c(paste(LETTERS[t1],collapse=''),paste(LETTERS[t2],collapse=''), paste(LETTERS[t3],collapse='')) } ans=unique(t(apply(pm,1,f))) print(head(ans,20),quote=FALSE) print(tail(ans,20),quote=FALSE) " 10人の子供(A,B,C,....,J)を3チームに割り当てる。 各チームは1人以上の子供が割り当てられる。 何通りの分け方があるか? 例. (A) (J) (B,C,D,E,F,G,H,I) (A,C) (H,I,J) (B,D,E,F,G) " rm(list=ls()) gr=expand.grid(1:8,1:8) gr1=gr[rowSums(gr)<10,] f=\(x) 10-sum(x) gr2=cbind(gr1,apply(gr1,1,f)) gr3=unique(t(apply(gr2,1,sort))) ans=apply(gr3,1,\(x) choose(10,x[1])*choose(10-x[1],x[2])/ length(unique(x)) ) sum(ans) res=cbind(gr3,ans) colnames(res)=NULL rownames(res)=NULL res 多摩南部地域病院に入院した事あるけど看護師がクソ過ぎ 患者のことはどうでもよく自分達の仕事さえまわせればそれでいいって感じ 底意地の悪いクソ女看護師の清○桃○は要注意 そのうえ入院環境もクソ過ぎ 狭いトイレ、備品等もボロいものばかりで配慮が足りなさすぎる 病院自体のシステムも不親切だし二度と行きたくない こんな病院は無いほうがいい 開業医もこんな病院紹介しないでほしい ganrikinto<-function(A=NULL,I=NULL,N=NULL,x=NULL){ # A:Principal I:Annual Interest N:months of repayment x:monthly repayment if(is.null(I)){ I=uniroot(\(m) x-A*m/(1-(1+m)^(-N)),c(1e-12,1-1e-12))$root*12 return(c('Annual Interest'=I)) } m=I/12 if(is.null(x)){ x=A*m/(1-(1+m)^(-N)) return(c('monthly repayment'=round(x,2))) } if(is.null(N)){ N=ifelse(x>A*m,log(x/(x-A*m))/log(m+1),Inf) return(c('repayment period in months'=round(N,2))) } if(is.null(A)){ A=x/m*(1-(1+m)^(-N)) return(c('Principal'=round(A))) } } f=\(n){ n1=n%%10 i=1 flg <- (n1*n1^i)%%10==n1 if(flg){ p=i }else{ while(!flg & i<2024){ i=i+1 flg <- (n1*n1^i)%%10==n1 } p=i } j=ifelse(2024%%p,2024%%p,p) (n1^p)%%10 } sapply(1:2024,f) |> table() 2024^2024 π^2024 の先頭の数字を算出 臨床医の嗜み:R言語での一般解を算出するコード calc=function(m,n) floor(10^(n*log10(m) -floor(n*log10(m)))) calc(2024,2024) calc(pi,2024) " 地球の公転周期は秒単位だと365日5時間48分46秒であるという 5時間48分46秒=20926秒 閏年は400年に97回 1日=86400秒 " p=20926/86400 ; p 97/400 - p 8/33 - p fn=\(x){ xp=x*p fl=floor(xp) ce=ceiling(xp) dfl=abs(fl/x-p) dce=abs(ce/x-p) c(min(dfl,dce),ifelse(dfl<dce,fl,ce),x) } dat=t(sapply(1:40,fn)) rbind(dat[which.min(dat[,1]),]) dat=t(sapply(1:10000,fn)) rbind(dat[which.min(dat[,1]),]) n=seq(4,400,4) dat=t(sapply(n,fn)) rbind(dat[which.min(dat[,1]),]) " 短径20cm 長径24cmの楕円のピザを楕円弧の長さが等しいように3つの楕円扇形に分割する。 楕円扇形の面積が最大になるようにするにはどのように切ればよいか? そのときの面積はいくつか?(答は整数値でよい)。 あらゆるリソースを用いてよい。 " a=12 b=10 n=3 S=\(t1,t2) a*b/2*(t2-t1) L=\(t1,t2) integrate(\(x) sqrt((-a*sin(x))^2+(b*cos(x))^2),t1,t2,abs.tol = 1e-24)$value Lt=L(-pi,pi) ; Lt # 全周長 t12t2=\(t1) uniroot(\(x) L(t1,x)-Lt/n,c(t1,t1+pi),tol=1e-24)$root t12S=\(t1) S(t1,t12t2(t1)) t12S=Vectorize(t12S) curve(t12S,0,pi) opt=optimize(t12S,c(-pi/2,pi/2),maximum = TRUE,tol=1e-24) ; opt source('toolmini.R') Plot(-a,a,zero=FALSE) Ellip(0i,a,b,lty=3) t1=opt$maximum ; t2=t12t2(t1) tt=seq(t1,t2,length=100) lines(a*cos(tt),b*sin(tt),lwd=2) T1=a*cos(t1)+1i*b*sin(t1) T2=a*cos(t2)+1i*b*sin(t2) seg(0i,T1,lty=2,lwd=2) seg(0i,T2,lty=2,lwd=2) Angle(T1,0i,T2) (pi*a*b - opt$obj)/2 S(t2,pi) rm(list=ls()) par(bty='l') calc=function(m=2024,n=2024) floor(10^(n*log10(m) -floor(n*log10(m)))) calc(2024,2024) calc(pi,2024) calc(exp(1),2024) " nを1以上2024以下の整数とする。 n^2024の先頭に現れる数字はなにが一番多いか? 頻度順に並べよ。 " calc.top =\(n) sapply(1:n,calc) |> table() |> sort() |> rev() re=calc.top(2024) re p=re/sum(re) calc.top(1000) sapply(1:1000,\(x) as.numeric(names(calc.top(x)[1]))) k=1:9 plot(k,p,pch=16) lines(k,log10(1+1/k),type='h',col=2) # ベンフォードの法則 医者板では全く相手にされないから数学板で高校生相手にイキれると思ってる尿瓶チンパンジジイ滑稽だね これが現実 1次方程式もできないド底辺特殊シリツ医大卒の記録 http://imagizer.imag...g923/2715/RosCsf.jpg 何度読んでも馬鹿すぎる。 男女別の割合と全体での割合から男女比が計算できるとも思わないとは。 なんでこんなのが大学に入れるわけよ? 裏口入学以外に説明がつく? 中学生でも解ける一次方程式の問題だろ。 それすらできない馬鹿が自信を持って発言。 >患者の男女比が必要なのもわからないのか? だとさ。 http://imagizer.imag...g923/9687/zNivZW.jpg read.cgi ver 07.5.4 2024/05/19 Walang Kapalit ★ | Donguri System Team 5ちゃんねる