臨床統計もおもしろいですよ、その3
1 :卵の名無しさん :2020/03/05(木) 20:17:05.03 ID:naSB8128.net 内科認定医受験の最低限の知識、 製薬会社の示してくる臨床データ、 論文の考察、 論文を書くときの正当性、 というのが、臨床統計の今までの目的の大きい部分でしたが、 AI=機械学習の基本も、結局は統計学と確率に支配されます。 そういう雑多な話をするスレです。 ※前スレ 臨床統計もおもしろいですよ、その1 https://egg.2ch.net/test/read.cgi/hosp/1493809494/ 臨床統計もおもしろいですよ、その2 https://egg.5ch.net/test/read.cgi/hosp/1540905566/
2 :卵の名無しさん :2020/03/05(木) 20:17:46.76 ID:naSB8128.net 100人の集団(クラスター)で感染者は1人として1日に延べ10回・人と接触し、 1人1回あたりの感染確率を1%、感染期間30日、潜伏期5日として SEIRモデルで計算すると、感染のピークは110日めでとても1〜2週間が山場とは言えない。 > SEIR2(contact_rate=10,transmission_probability=0.01 + ,infectious_period=30,latent_period=5,mu=0, + nu=0, s=99,e=0,i=1,r=0,timepoints = seq(0,365,by=0.5),axes=TRUE) Ro = 3 peak time I = 109.5 peak time E = 89 グラフにすると https://i.imgur.com/EiG2HzI.jpg
3 :卵の名無しさん :2020/03/05(木) 22:34:57.14 ID:vGY9zXsV.net https://i.imgur.com/obIkECO.jpg
4 :卵の名無しさん :2020/03/06(金) 09:29:34.46 ID:H1NcHDj6.net 病的な虚言癖と妄想癖の精神科医 古根高の病名を診断するスレ https://egg.5ch.net/test/read.cgi/hosp/1529634250/ 古根高(サッポロファクトリーメンタルクリニック院長)(五浪で底辺医大に進学した知恵遅れ)のご尊顔 https://i.imgur.com/u8rBskz.jpg トラブルメーカーゆえ市立札幌病院を解雇された過去あり https://i.imgur.com/upWjhiA.jpg
5 :卵の名無しさん :2020/03/07(土) 09:19:47.77 ID:jrmaVufM.net 週末はこれで遊べる。 nCov2019 for studying COVID-19 coronavirus outbreak https://guangchuangyu.github.io/nCov2019/
6 :卵の名無しさん :2020/03/07(土) 10:32:32.16 ID:jrmaVufM.net >>5 まず、Rを最新版に更新して if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("BiocStyle") remotes::install_github("GuangchuangYu/nCov2019", dependencies = TRUE) が必要だった。
7 :卵の名無しさん :2020/03/07(土) 10:47:18.70 ID:jrmaVufM.net BiocStyle-defunct {BiocStyle} R Documentation Defunct functions in package ‘BiocStyle’ Description These functions are defunct and no longer available. Details The following functions are no longer available; use the replacement indicated below: latex_old, latex2: latex pdf_document_old, pdf_document2: pdf_document html_document_old, html_document2: html_document
8 :卵の名無しさん :2020/03/07(土) 15:19:10.32 ID:KnbHZezX.net >>1 受験を控えている高校生です 他のスレでアナタはコンプ性Fラン症候群という不治の病であるとお聞きしました 高学歴の私と違ってFランであるアナタは今はどのようなお仕事をされてるんですか?
9 :卵の名無しさん :2020/03/07(土) 19:47:32.16 ID:jrmaVufM.net Coronavirus COVID-19 outbreak statistics and forecast http://www.bcloud.org/e/ これは使えるな http://www.bcloud.org/e/session/a9b2a4a3e4ef04951a3e45743536a8f3/download/dataDownloadWorld?w=
10 :卵の名無しさん :2020/03/08(日) 14:49:32.83 ID:WlKU3B2d.net SEIRモデルで有病率を1%に固定して、集団のサイズを変化させてシミュレーションしてみたけどピークは変わらないな。 このモデルでは集会規模の大小には影響されないということになるな。 https://i.imgur.com/343K91V.png 有病率を変化させて流行の変遷をグラフにすると、 https://i.imgur.com/SZ15LKT.png https://i.imgur.com/gnJVFnd.png 有病率を40%くらいに引き上げるとオリンピックのときには流行が収束していることになるwwww
11 :卵の名無しさん :2020/03/08(日) 19:19:22.20 ID:WlKU3B2d.net rm(list=ls()) graphics.off() f <- function(n,p) 1-(1-p)^n vf=Vectorize(f) n=seq(1,1000,by=1) plot(n,vf(n,0.001),bty='l') abline(h=0.5,lty=3) " 藤井聡が 集会の参加人数と感染拡大確率の表を公開している。 ://i.imgur.com/f3Tes5g.jpg 同じアルゴリズムを使って有病率が1/10000のときに感染拡大確率を0.05未満にしたい。 何人以上の集会を禁じることによってそれが達成できるか計算せよ。 " F <- function(p=0.001,n=1000,p0=0.5){ # p=有病率,n:最大参加人数,p0:感染拡大確率 sub <- function(n,p) 1-(1-p)^n uniroot(function(x) sub(x,p)-p0,c(0,n))$root } F() F(p=1e-4,n=1e4,p0=0.05)
12 :卵の名無しさん :2020/03/08(日) 19:19:38.25 ID:WlKU3B2d.net " p0= 1-(1-p)^n (1-p)^n=1-p0 n*log(1-p)=log(1-p0) n=log(1-p0)/log(1-p) " pp02n <- function(p,p0) log(1-p0)/log(1-p) # p=有病率,p0:感染拡大確率 pp02n(1/1000,0.05) p=1/10^seq(1,5,by=0.01) p=rev(p) plot(p,pp02n(p,0.05),log='x',type='l',bty='l',ylab='参加者数',xlab='有病率', main='感染拡大確率別許容参加人数') lines(p,pp02n(p,0.10),lty=2,lwd=2) lines(p,pp02n(p,0.25),lty=3,lwd=3) lines(p,pp02n(p,0.50),lty=4,lwd=4) legend('top',bty='n',legend=c(0.05,0.10,0.25,0.5),lty=1:4,lwd=1:4) " p0= 1-(1-p)^n (1-p)^n=1-p0 log(1-p)=log(1-p0)/n 1-p=exp(log(1-p0)/n) p=1-exp(log(1-p0)/n) " n=465 p0=0.05 np02p <- function(n,p0) 1-exp(log(1-p0)/n) np02p(465,0.05)
13 :卵の名無しさん :2020/03/08(日) 20:01:02.96 ID:WlKU3B2d.net " 東京都は21日、新型コロナウイルスの感染拡大を防ぐため、22日から3月15日までの3週間、都が主催する500人以上の大規模な屋内イベントは、原則延期か中止にする方針を発表した. " by=1e-5 p0=seq(by,0.5,by) plot(np02p(500,p0),p0,log='x',bty='l',type='l',lwd=2,main='500人参加集会', xlab='有病率(log)',ylab='感染拡大確率')
14 :卵の名無しさん :2020/03/09(月) 08:18:17.40 ID:cZTNWJLE.net Same Stats, Different Graphs: Generating Datasets with Varied Appearance and Identical Statistics through Simulated Annealing Justin Matejka and George Fitzmaurice Autodesk Research, Toronto Ontario Canada {first.last}@autodesk.com
15 :卵の名無しさん :2020/03/09(月) 08:21:02.64 ID:cZTNWJLE.net https://blogs.sas.com/content/iml/2019/04/17/create-version-of-anscombes-quartet.html
16 :卵の名無しさん :2020/03/09(月) 08:23:47.32 ID:cZTNWJLE.net >>14 https://www.researchgate.net/profile/Justin_Matejka/publication/ 316652618_Same_Stats_Different_Graphs_Generating_Datasets_with_Varied_Appearance_and_Identical_Statistics_through_Simulated_Annealing/links/ 599aef990f7e9b892bacff30/Same-Stats-Different-Graphs-Generating-Datasets-with-Varied-Appearance-and-Identical-Statistics-through-Simulated-Annealing.pdf
17 :卵の名無しさん :2020/03/09(月) 08:24:55.87 ID:cZTNWJLE.net https://youtu.be/It4UA75z_KQ
18 :卵の名無しさん :2020/03/09(月) 19:01:23.07 ID:Ost2Ton5.net Sys.setenv(lang='en') rm(list=ls()) library(gtools) n=4 a=permutations(n,n) # nの順列 r=nrow(a) f<-function(x){ # x=c(1,2,3) -> rbind(a[1],a[2],a[3]) n=length(x) ans=NULL for(i in 1:n){ ans=rbind(ans,a[i,]) } return(ans) } (b=combn(r,n,f)) # すべての行列の配列 b[,,i] (c=dim(b)[3]) # その個数 diag2 <- function(x){ # 右上から左下への対角線の配列を返す n=nrow(x) ans=numeric(n) for(i in 1:n) ans[i]=x[i,n+1-i] return(ans) }
19 :卵の名無しさん :2020/03/09(月) 19:01:30.17 ID:Ost2Ton5.net g <- function(x){ # 列、対角線の要素がすべて異なればTRUEを返す n=nrow(x) if(length(unique(diag (x))) < n) return(FALSE) if(length(unique(diag2(x))) < n) return(FALSE) flag=TRUE for(i in 1:n){ if(length(unique(x[,i])) < n){ flag=FALSE break } } return(flag) } count=0 for(i in 1:c){ count=count+g(b[,,i]) } count
20 :卵の名無しさん :2020/03/09(月) 20:15:13.42 ID:Ost2Ton5.net >>18 これバグがあるな
21 :卵の名無しさん :2020/03/09(月) 21:27:07.44 ID:Ost2Ton5.net >>20 バグ修正版 https://rio2016.5ch.net/test/read.cgi/math/1501755792/789
22 :卵の名無しさん :2020/03/10(火) 17:40:25.16 ID:Fd2PNckm.net n=1e3 p0=runif(n,0,1) oz0=p0/(1-p0) pLR=0.7/(1-0.9) #TP/FP nLR=0.3/0.9 # FN/TN oz1=oz0*nLR^2*pLR^2 p1=oz1/(1+oz1) quantile(p1,c(.025,0.5,.95))
23 :卵の名無しさん :2020/03/10(火) 18:13:44.48 ID:Fd2PNckm.net sn=0.7 sp=0.9 n=1e7 p0=runif(n,0,1) oz0=p0/(1-p0) pLR=sn/(1-sp) #TP/FP nLR=(1-sn)/sp # FN/TN oz1=oz0*nLR^2*pLR^2 p1=oz1/(1+oz1) BEST::plotPost(p1,showMode = T) HDInterval::hdi(p1) quantile(p1,c(.025,0.5,.975)) summary(p1) MAP <- function(x) { dens <- density(x) mode_i <- which.max(dens$y) mode_y <- dens$y[mode_i] c(x=mode_x, y=mode_y) } MAP(p1)
24 :卵の名無しさん :2020/03/11(水) 07:17:27.96 ID:w+Rtdm4+.net PCR2prob <- function( sn=0.7, # sensitivity sp=0.9, # specificity plus=2, # how many positive result? minus=2, # how many negative result? n=1e7, print=TRUE) # how large the simulation { p0=runif(n,0,1) oz0=p0/(1-p0) # prob -> odds pLR=sn/(1-sp) # TP/FP nLR=(1-sn)/sp # FN/TN oz1=oz0*pLR^plus*nLR^minus # Bayesian formula p1=oz1/(1+oz1) # odds -> prob if(print){ BEST::plotPost(p1,showMode =T) # show mode instead of mean print(HDInterval::hdi(p1)) # Highest Density Interval print(quantile(p1,c(.025,0.5,.975))) # 95%CI by quantile print(summary(p1)) # mean, median } invisible(p1) } PCR2prob(n=1e5,minus=3) PCR2prob(n=1e5,minus=4)
25 :卵の名無しさん :2020/03/11(水) 07:17:32.52 ID:w+Rtdm4+.net MAP <- function(x) { dens <- density(x) mode_i <- which.max(dens$y) mode_x <- dens$x[mode_i] mode_y <- dens$y[mode_i] c(x=mode_x, y=mode_y) } MAP(p1)['x'] # show mode f <- function(x) mean(PCR2prob(minus=x,n=1e4,print=F)) plot(sapply(1:10,f),bty='l',pch=19) abline(h=c(0.05,0.1,0.5),col='gray',lty=3)
26 :卵の名無しさん :2020/03/11(水) 10:07:10.64 ID:NmX2BxF2.net PCR2prob <- function( sn=0.7, # sensitivity sp=0.9, # specificity plus=2, # how many positive result? minus=2, # how many negative result? n=1e7, # how large the simulation p0=runif(n,0,1), print=FALSE) { oz0=p0/(1-p0) # prob -> odds pLR=sn/(1-sp) # TP/FP nLR=(1-sn)/sp # FN/TN oz1=oz0*pLR^plus*nLR^minus # Bayesian formula p1=oz1/(1+oz1) # odds -> prob if(print){ BEST::plotPost(p1,showMode =T) # show mode instead of mean print(HDInterval::hdi(p1)) # Highest Density Interval print(quantile(p1,c(.025,0.5,.975))) # 95%CI by quantile print(summary(p1)) # mean, median } invisible(p1) } mean(PCR2prob(n=1e5,minus=3)) PCR2prob(n=1e5,minus=4)
27 :卵の名無しさん :2020/03/11(水) 10:55:05.81 ID:NmX2BxF2.net PCR2prob <- function( sn=0.7, # sensitivity sp=0.9, # specificity plus=2, # how many positive result? minus=2, # how many negative result? n=1e7, # how large the simulation p0=rbeta(n,1,1), # prior distribution print=TRUE) { oz0=p0/(1-p0) # prob -> odds pLR=sn/(1-sp) # TP/FP nLR=(1-sn)/sp # FN/TN oz1=oz0*pLR^plus*nLR^minus # Bayesian formula p1=oz1/(1+oz1) # odds -> prob if(print){ BEST::plotPost(p1,showMode =T) # show mode instead of mean print(HDInterval::hdi(p1)) # Highest Density Interval print(quantile(p1,c(.025,0.5,.975))) # 95%CI by quantile print(summary(p1)) # mean, median } invisible(p1) # return posterior distribution }
28 :卵の名無しさん :2020/03/11(水) 15:05:50.80 ID:w+Rtdm4+.net sim <- function(m=2){ count=0 i=0 while(count<m){ i=i+1 count = count + sample(6,1)==1 } return(i) } k=1e5 re=replicate(k,sim()) tbl=table(re) ; tbl which.max(tbl) plot(tbl/k,bty='l')
29 :卵の名無しさん :2020/03/11(水) 17:01:09.37 ID:w+Rtdm4+.net "1が累計m(=2)回出るまでサイコロを振って、振った回数を当てるギャンブルがある。 何回目に賭けるのがベストか?" sim <- function(m=2){ pip1=0 # 1の目の出た回数 i=0 # サイコロを振った回数 while(pip1 < m){ i=i+1 pip1 = pip1 + (sample(6,1)==1) } return(i) } k=1e5 re=replicate(k,sim(3)) tbl=table(re) ; tbl which.max(tbl) plot(tbl/k,bty='l')
30 :卵の名無しさん :2020/03/11(水) 17:01:26.29 ID:w+Rtdm4+.net # bg <- function(x,print=FALSE){ # big gambling f <- function(n,m=x,p=1/6) choose(n-1,m-1)*p^(m-1)*(1-p)^(n-m)*p nn=1:(10*x) y=optimize(function(n) f(n),nn,maximum=TRUE)$maximum if(print){ plot(nn,sapply(nn,f),bty='l',pch=19) yy=c(floor(y),ceiling(y)) cat(c(f(yy[1]),f(yy[2])),'\n') } return(floor(y)) } x=2:100 y=sapply(x,bg) plot(x,y,bty='l') (lm=lm(y~x)) summary(lm) abline(lm) y <- function(x) 6*x-6 y(1000)
31 :卵の名無しさん :2020/03/12(木) 20:47:51.16 ID:Hoqw8S6P.net Rに移植 allocate.rooms <- function(m,n){ # m:rooms n:people if(m==n) return(factorial(m)) else if(m==1) return(1) else m*Recall(m,n-1) + m*Recall(m-1,n-1) } Cに移植 #include<stdio.h> long factorial(long n) { long re = 1; long k; for(k=1;k <=n;k++) {re *= k;} return re; } long rooms(int m, int n){ if(m==n) { return factorial(m);} else if(m==1){ return 1;} else{ return m * rooms(m,n-1) + m * rooms(m-1,n-1); } } void main( int argc, char *argv[] ){ int m,n; long ways; m=atoi(argv[1]); n=atoi(argv[2]); ways=rooms(m,n); printf("%d\n",ways); }
32 :卵の名無しさん :2020/03/12(木) 20:48:19.73 ID:Hoqw8S6P.net m部屋n人空部屋なしの場合の数で a[m,n] = a[m,n-1]×m + a[m-1,n-1]×m *Main> let s m n = if m == n then (product [1..m]) else if m==1 then 1 else m*(s (m) (n-1)) + m*(s (m-1) (n-1)) *Main> s 6 12 953029440
33 :卵の名無しさん :2020/03/14(土) 12:34:11.45 ID:E6fZCvx4.net ド底辺シリツ医大の新入生100人のうち10人は学力考査で合格した正規入学生で残りは裏口入学生であるとする。 誰が正規入学かを確かめる必要がある。審議官の女医(嘘つきかどうかは不明)が無作為に新入生を呼び出して 「あなたのいうことが正しければ手コキかフェラをしてあげる、そうでなければセンズリを命じる」という試問をする。 論理的思考ができる正規入学生はこの質問に回答してフェラをしてもらっているが、裏口入学生は正解がだせずにセンズリを命じられている。 正規入学生10人を同定できたら一連の試問は終了とする。 センズリを命じられる裏口入学生の期待値を求めよ。 rm(list=ls()) students=rep(1:2,c(10,90)) picked=NULL flag=FALSE sim <- function(){ while(flag==FALSE){ i=sample(length(students),1) picked=c(picked,students[i]) students=students[-i] flag=sum(picked==1)==10 } sum(picked==2) } k=1e5 mean(replicate(k,sim())) 900/11
34 :卵の名無しさん :2020/03/14(土) 14:06:57.29 ID:E6fZCvx4.net " 検査陽性 検査陰性 カレー頻食 a(=18) b(=30) カレー稀食 c(=32) d(=20) Usage riskratio(X, Y, m1, m2, conf.level=0.95, p.calc.by.independence=TRUE) Arguments X The number of disease occurence among exposed cohort. # a Y The number of disease occurence among non-exposed cohort. # c m1 The number of individuals in exposed cohort group. # a+b m2 The number of individuals in non-exposed cohort group. # c+d " curry <- function( prev=0.10, # prevalence N=1, # population sn=0.7, # sensitivity sp=0.9, # specificity a=18,b=30, c=32,d=20) { PPV=N*prev*sn/(N*prev*sn+N*(1-prev)*(1-sp)) NPV=N*(1-prev)*sp/(N*(1-prev)*sp+N*prev*(1-sn)) library(fmsb) p0=riskratio(a,c,a+b,c+d)$p.value p1=riskratio(a*PPV,c*PPV,a*PPV+b*NPV,c*PPV+d*NPV)$p.value c(p0,p1) } curry()
35 :卵の名無しさん :2020/03/14(土) 14:27:13.93 ID:E6fZCvx4.net 某国の新型コロナ感染症の有病率を0.1、PCR検査の感度を0.7 特異度を0.9とする 検査陽性陰性の人を無作為に50人ずつ集めて、カレーを頻回に食べているかを調査した結果が 以下の通りであった。 検査陽性 検査陰性 カレー頻食 a(=18) b(=30) カレー稀食 c(=32) d(=20) カレーを頻食すると新型コロナ感染に罹りにくいと結論できるか? 危険率0.05で検定せよ。
36 :卵の名無しさん :2020/03/14(土) 14:33:35.06 ID:E6fZCvx4.net >>35 有病率によって有意かどうか判断が変わるのは興味ぶかい。 https://i.imgur.com/qjjIjnm.png
37 :卵の名無しさん :2020/03/14(土) 20:12:22.41 ID:E6fZCvx4.net 症数例から断定的な判断をするのは危険なので 周辺度数を固定してMCMCしてp値の信頼区間を出してみた。 # prev=0.10 # prevalence sn=0.7 # sensitivity sp=0.9 # specificity a=18 b=30 c=32 d=20 (x=matrix(c(a,b,c,d),2, byrow=TRUE)) N=sum(x) PPV=N*prev*sn/(N*prev*sn+N*(1-prev)*(1-sp)) NPV=N*(1-prev)*sp/(N*(1-prev)*sp+N*prev*(1-sn)) library(BayesFactor) bf=contingencyTableBF(x,sampleType = 'indepMulti',fixedMargin = 'cols', posterior = TRUE, iteration=1e4) head(bf) N=sum(x) a=N*bf[,1] b=N*bf[,2] c=N*bf[,3] d=N*bf[,4] library(fmsb) p.value=riskratio(a*PPV,c*PPV,a*PPV+b*NPV,c*PPV+d*NPV)$p.value BEST::plotPost(p.value,compVal = 0.05) https://i.imgur.com/w2RzDLB.png p値の期待値は0.152 有意水準0.05で有意差ありとされる確率は約4割という結果がえられた。
38 :卵の名無しさん :2020/03/14(土) 20:16:36.23 ID:E6fZCvx4.net >>37 タイプミス修正 症数例から ↓ 少数例
39 :卵の名無しさん :2020/03/15(日) 08:45:40.74 ID:LJXDSnIu.net https://id.fnshr.info/2019/07/20/r3-6-0/ で知った 例 rep("1:6",5) paste(rep("1:6",5),collapse=',') str2lang(paste("expand.grid(",paste(rep("1:6",5),collapse=','),')')) eval(str2lang(paste("expand.grid(",paste(rep("1:6",5),collapse=','),')'))) (s=rep("1:6",5)) (str=paste(s,collapse=',')) (lang=str2lang(paste("expand.grid(",str,')'))) eval(lang)
40 :卵の名無しさん :2020/03/15(日) 08:53:22.23 ID:LJXDSnIu.net fn <− function(n){ library(gmp) gr=eval(str2lang(paste(Wexpand.grid(W,paste(rep(W1:6W,n),collapse=V,V),W)W))) r=as.bigq(sum(apply(gr,1,function(x) prod.bigz(x)%%6==0)))/as.bigq(nrow(gr)) r2=capture.output(r)[2] substr(r2,5,nchar(r2)) } for(i in 1:9) cat(i,V:V,fn(i),V\nV) # str2lang str2expression ?str2lang rep(W1:6W,5) paste(rep(W1:6W,5),collapse=V,V) str2lang(paste(Wexpand.grid(W,paste(rep(W1:6W,5),collapse=V,V),V)V)) eval(str2lang(paste(Wexpand.grid(W,paste(rep(W1:6W,5),collapse=V,V),V)V))) (s=rep(W1:6W,5)) (str=paste(s,collapse=V,V)) (lang=str2lang(paste(Wexpand.grid(W,str,V)V))) eval(lang)
41 :卵の名無しさん :2020/03/15(日) 16:22:13.79 ID:5ItQZoKJ.net >>1 受験を控えている女子高校生です 他のスレでアナタはコンプ性Fラン症候群という不治の病であるとお聞きしました 高学歴の私と違ってFランであるアナタは今はどのようなお仕事をされてるんですか?
42 :卵の名無しさん :2020/03/15(日) 18:08:53.16 ID:LJXDSnIu.net >41 高学歴ならこれに答えてみ! 2020/03/13 11:00時点で東京都で1524人検査して87人陽性と報告されている # https://data-science.gr.jp/theory/tpd_negative_binomial_distribution.html 発見率を87/1524は一定仮定して 新型コロナ陽性患者を10人集めたいとする。 必要な被検者数の期待値と95%信頼区間を求めよ。 # Monte Carlo p=87/1524 failure=0:500 success=10 fail=rnbinom(1e6,success,p) hist(fail,freq=F,col='seashell2',main='',xlab='failures',ylab='probability',breaks=30) lines(failure,dnbinom(failure,success,p),lwd=2) mean(fail) + success success*(1-p)/p + success # mean(fail)=kq/p HDInterval::hdi(x) quantile(x,c(0.025,0.5,0.975)) qnbinom(c(0.025,0.975),success,p)
43 :卵の名無しさん :2020/03/15(日) 18:09:31.68 ID:LJXDSnIu.net シミュレーション解の方が回数上限なしでプログラムできるな。 # Simulation sim <- function(){ i=0 s10=0 while(s10!=10){ i=i+1 s10 = s10 + sample(1:0,1,prob=c(p,1-p)) s10 == 10 } i } k=1e4 re=replicate(k,sim()) mean(re) HDInterval::hdi(re)
44 :卵の名無しさん :2020/03/15(日) 18:33:54.47 ID:LJXDSnIu.net 2020/03/13 11:00時点で東京都で1524人検査して87人陽性と報告されている # https://data-science.gr.jp/theory/tpd_negative_binomial_distribution.html 発見率は平均値87/1524、0.01-0.10の区間の二項分布に従うと仮定する。 新型コロナ陽性患者を10人集めたい。 必要な被検者数の期待値と95%信頼区間を求めよ。 library(BayesFactor) bf=BayesFactor::proportionBF(87,1524, p, posterior = TRUE,iter=1e4,nullInterval = c(0.01,0.10)) head(bf) pp=(bf[,'p']) plot(pp,bty='l') mean(pp) BEST::plotPost(pp,showMode = T) qq=1-p r=10 trials=r*qq/pp + r BEST::plotPost(trials) HDInterval::hdi(trials)
45 :卵の名無しさん :2020/03/15(日) 18:40:27.27 ID:LJXDSnIu.net >>41 ド底辺頭脳でも数くらい数えられるんだろ? 1つのサイコロを10000回投げ、出たの目の積が6の倍数になる確率を分数でもとめよ。
46 :卵の名無しさん :2020/03/15(日) 19:19:43.78 ID:5ItQZoKJ.net >>45 どうでもいいから早くコンプ薬屋を呼んで来なよ そんなだからモテないんだお by 都内高学歴JK
47 :卵の名無しさん :2020/03/15(日) 19:31:13.08 ID:LJXDSnIu.net >>46 数も数えられんの?ガイジ?
48 :卵の名無しさん :2020/03/15(日) 21:05:46.74 ID:LJXDSnIu.net >>35 curry <- function( prev=0.10, # prevalence sn=0.7, # sensitivity sp=0.9, # specificity a=18,b=30, c=32,d=20) { N=1 # population PPV=N*prev*sn/(N*prev*sn+N*(1-prev)*(1-sp)) NPV=N*(1-prev)*sp/(N*(1-prev)*sp+N*prev*(1-sn)) library(fmsb) p0=riskratio(a,c,a+b,c+d)$p.value a=a*PPV+b*(1-NPV) b=b*NPV+a*(1-PPV) c=c*PPV+d*(1-NPV) d=d*NPV+c*(1-PPV) p1=riskratio(a,c,a+b,c+d)$p.value c(p0,p1) } curry() prevs=seq(0.01,0.99,by=0.01) plot(prevs,sapply(prevs,function(x) curry(prev=x)[2]),bty='l',type='l',lwd=2, ylim=c(0,1),ylab='p-value',xlab='prevalence') curry(a=36,c=64,b=60,d=40)
49 :卵の名無しさん :2020/03/15(日) 21:06:08.12 ID:LJXDSnIu.net prev=0.10 # prevalence sn=0.7 # sensitivity sp=0.9 # specificity a=18 b=30 c=32 d=20 N=sum(x) PPV=N*prev*sn/(N*prev*sn+N*(1-prev)*(1-sp)) NPV=N*(1-prev)*sp/(N*(1-prev)*sp+N*prev*(1-sn)) a=a*PPV+b*(1-NPV) b=b*NPV+a*(1-PPV) c=c*PPV+d*(1-NPV) d=d*NPV+c*(1-PPV) (x=matrix(c(a,b,c,d),2, byrow=TRUE)) library(BayesFactor) bf=contingencyTableBF(x,sampleType = 'indepMulti',fixedMargin = 'cols', posterior = TRUE, iteration=1e4) head(bf) N=sum(x) a=N*bf[,1] b=N*bf[,2] c=N*bf[,3] d=N*bf[,4] library(fmsb) a=a*PPV+b*(1-NPV) b=b*NPV+a*(1-PPV) c=c*PPV+d*(1-NPV) d=d*NPV+c*(1-PPV) p=riskratio(a,c,a+b,c+d)$p.value p.value=as.vector(p) BEST::plotPost(p.value,compVal = 0.05)
50 :卵の名無しさん :2020/03/15(日) 21:10:18.08 ID:LJXDSnIu.net F-value = 2/(1/sensitivity + 1/PPV ) 感度とPPVの調和平均
51 :卵の名無しさん :2020/03/16(月) 00:14:20.50 ID:RRM1UH+L.net # n個からr個を選んで得られる順列の総数をP(n, r)とする. 任意のr>1に対して, P(n, r)は平方数でないことを示せ. rm(list=ls()) library(gmp) library(Rmpfr) fn <- function(n){ r=2:n a=chooseZ(n,r)*factorialZ(r) # nPr b=mpfr(a,1e5) re=is.whole(sqrt(b)) r[which(re)] } fn(1000)
52 :卵の名無しさん :2020/03/16(月) 06:35:35.55 ID:7zdTC761.net >>51 早くFラン仲間のコンプ薬屋を呼んで来なよ お勉強してくれるって言ってるお by 都内高学歴JK
53 :卵の名無しさん :2020/03/16(月) 06:40:26.18 ID:RRM1UH+L.net >>52 ド底辺頭脳でも数くらい数えられるんだろ? 1つのサイコロを10000回投げ、出たの目の積が6の倍数になる確率を分数でもとめよ。
54 :卵の名無しさん :2020/03/16(月) 07:32:37.13 ID:RRM1UH+L.net 国内で患者数が大幅に増えたときに備えた医療提供体制の確保について 今後、国内で新型コロナウイルス感染症患者数が大幅に増えたときに備え、各都道府県、保健所設置市及び特別区(以下「都道府県等」という。)における外来を受診する患者数等について、 以下の数式を用いて計算いただき、ピーク時の医療需要の目安としてご活用の上、必要な医療提供体制を確保していただくようお願いいたします。 (1)(ピーク時において1 日あたり新たに新型コロナウイルス感染症を疑って外来を受診する患者数)= (0-14 歳人口)×0.18/100+(15-64 歳人口) ×0.29/100+(65 歳以上人口) ×0.51/100 (2)(ピーク時において1 日あたり新型コロナウイルス感染症で入院治療が必要な患者数)= (0-14 歳人口)×0.05/100+(15-64 歳人口)×0.02/ 100+(65 歳以上人口) ×0.56/100 (3)(ピーク時において1 日あたり新型コロナウイルス感染症で重症者として治療が必要な患者数)= (0-14 歳人口)×0.002/100+(15-64 歳人口) ×0.001/100+(65 歳以上人口) ×0.018/100 注1)ピーク時は、各都道府県等において疫学的関連性が把握できない程度に感染が拡大した時点から概ね3か月後に到来すると推計されている。ただし、公衆衛生上の対策を行うことにより、ピークが下がるとともに後ろ倒しされる。 注2)重症者とは、集中治療や人工呼吸器を要する管理が必要な患者を指す。 注3)当該計算式は、都道府県等の単位以下における医療提供体制を確保するためのものであるとともに、各都道府県等によってピークを迎える時期が異なるため、 全国の人口を用いて計算することや単純に各自治体が算出するピークの数値を足し合わせることは、不適切な取扱いとなることに留意いただきたい。なお、当該計算式については、 今後新たな知見等により変更される可能性がある。 注4)実際には、ピーク時に至るまでの日々の患者数の増加はばらつきがあり、増加曲線は推計通りの形にならない可能性が高いため、 現実の患者の発生動向も踏まえて適切に体制を確保することが必要。 注5)当該計算式については、今後新たな知見等により変更される可能性がある。
55 :卵の名無しさん :2020/03/16(月) 07:32:47.99 ID:RRM1UH+L.net v=c(0.18,0.29,0.51, 0.05,0.02,0.56, 0.002,0.001,0.018) (mat=matrix(v*100,3,byrow=T))
56 :卵の名無しさん :2020/03/16(月) 07:35:19.73 ID:7zdTC761.net 早くFラン仲間のコンプ薬屋を呼んで来なよ 童貞仲間でしょぉお by 都内高学歴JK
57 :卵の名無しさん :2020/03/16(月) 07:54:48.43 ID:RRM1UH+L.net v=c(0.18,0.29,0.51, 0.05,0.02,0.56, 0.002,0.001,0.018) (mat=matrix(v*100,3,byrow=T)) # https://www.toukei.metro.tokyo.lg.jp/juukiy/2019/jy19qa0200.xls (x=matrix(c(1601348,9035668,3103714)/1e4,ncol=1)) round(mat%*%x) > round(mat%*%x) [,1] [1,] 44915 [2,] 19989 [3,] 681 重症者の5割が死亡するとすると、都内で340人の死者予想。
58 :卵の名無しさん :2020/03/16(月) 07:56:56.15 ID:RRM1UH+L.net >>56 コンプはこれに答がだせなかったんだが、あんたは出せる? 解答できなきゃ、コンプと同レベルの頭脳と認定してあげよう。 24時間待ってあげる。 問題 某国の新型コロナ感染症の有病率を0.1、PCR検査の感度を0.7 特異度を0.9とする 検査陽性陰性の人を無作為に50人ずつ集めて、カレーを頻回に食べているかを調査した結果が 以下の通りであった。 検査陽性 検査陰性 カレー頻食 a(=18) b(=30) カレー稀食 c(=32) d(=20) カレーを頻食すると新型コロナ感染に罹りにくいと結論できるか? 危険率0.05で検定せよ。
59 :卵の名無しさん :2020/03/16(月) 08:12:57.28 ID:RRM1UH+L.net 国内で患者数が大幅に増えたときに備えた医療提供体制の確保について 今後、国内で新型コロナウイルス感染症患者数が大幅に増えたときに備え、各都道府県、保健所設置市及び特別区(以下「都道府県等」という。)における外来を受診する患者数等について、 以下の数式を用いて計算いただき、ピーク時の医療需要の目安としてご活用の上、必要な医療提供体制を確保していただくようお願いいたします。 (1)(ピーク時において1 日あたり新たに新型コロナウイルス感染症を疑って外来を受診する患者数)= (0-14 歳人口)×0.18/100+(15-64 歳人口) ×0.29/100+(65 歳以上人口) ×0.51/100 (2)(ピーク時において1 日あたり新型コロナウイルス感染症で入院治療が必要な患者数)= (0-14 歳人口)×0.05/100+(15-64 歳人口)×0.02/ 100+(65 歳以上人口) ×0.56/100 (3)(ピーク時において1 日あたり新型コロナウイルス感染症で重症者として治療が必要な患者数)= (0-14 歳人口)×0.002/100+(15-64 歳人口) ×0.001/100+(65 歳以上人口) ×0.018/100 問題 : ある都市でピーク時に外来1000人、入院600人、重症20人であったとすると、この都市の14歳以下の人口は何人と推測されるか?
60 :卵の名無しさん :2020/03/16(月) 13:47:58.64 ID:ulqLVmOp.net >>58 評価するには症例数が2桁少ないですわね それはさておき昔の教授選では積み上げたペーパーの高さを比べたと聞きますが Fランのアナタは採択されたペーパーの1本もあるのかしら?W by 都内Sラン女子高生
61 :卵の名無しさん :2020/03/16(月) 14:14:20.21 ID:ncInpF91.net やったー 国試無事合格!
62 :卵の名無しさん :2020/03/16(月) 23:11:32.64 ID:5fhQT99v.net >>60 足りないのは標本数でなくてオマエの頭だよ。
63 :卵の名無しさん :2020/03/16(月) 23:31:54.26 ID:5fhQT99v.net 数日前、科学技術省の共同予防および制御メカニズムの科学研究グループからの良いニュースがありました。 国家緊急予防および制御薬物工学技術研究センターと深セン第三人民病院LeレイおよびLi英Yチームによって完成された「ファピラビル」。 「新しいコロナウイルス肺炎(COVID-19)患者の安全性と有効性に関する臨床研究」は、予備的な臨床結果に達しています(登録番号: ChiCTR2000029600)。 研究により、ファビピラビルは、ウイルスクリアランスを促進することにより、新しいコロナウイルス肺炎の進行を緩和することが示されています。研究結果は、中国工学院のジ ャーナルに提出されました 合計80人の患者がこの臨床試験に登録され、35人の患者がファピラビルで治療され、対照群は年齢、性別、および疾患の重症度が治療群と一致した新冠動脈肺炎の患者45人 でした。 iタブレット治療。薬物投与からウイルス除去までの時間の中央値、治療の14日目の胸部画像の改善率、および安全性を2つのグループ間で比較しました。 結果は、ウイルス除去の観点から、治療後のファピラビル試験群の患者のウイルス除去時間の中央値(ウイルス核酸陰性)は、それぞれ4日と11日であった対照群のそれよりも 有意に短いことを示しました。別の重要な指標として、患者の胸部画像の改善、治療群と対照群の改善率はそれぞれ91.43%と62.22%でした。同時に、対照群と比較し て、フェラビル治療群は副作用が少なく、忍容性が良好でした。 2つの患者グループのベースライン特性はすべて同等でした。
64 :卵の名無しさん :2020/03/16(月) 23:35:45.65 ID:5fhQT99v.net >>63 PCRの方は中央値だけでデータがないので統計処理は辿れないが、比率の方が計算できる。 > round(35*0.9143) [1] 32 > round(45*0.6222) [1] 28 > > library(fmsb) > (mat=cbind(c(32,28),c(35,45)-c(32,28))) [,1] [,2] [1,] 32 3 [2,] 28 17 > chisq.test(mat) Pearson's Chi-squared test with Yates' continuity correction data: mat X-squared = 7.4667, df = 1, p-value = 0.006285 > fisher.test(mat) Fisher's Exact Test for Count Data data: mat p-value = 0.003669 alternative hypothesis: true odds ratio is not equal to 1 95 percent confidence interval: 1.59236 37.25281 sample estimates: odds ratio 6.33637
65 :卵の名無しさん :2020/03/16(月) 23:36:27.82 ID:5fhQT99v.net MCMCして信頼区間で比較すると、 https://i.imgur.com/3uZNIwI.png
66 :卵の名無しさん :2020/03/17(火) 00:11:08.71 ID:+SeXrvpf.net >>63 PCR陰転までの日数の中央値が アビガン群35例で4日、対照群で11日としても 次のような分布なら有効とは言い難い。 https://i.imgur.com/TjJshMu.png > summary(A) Min. 1st Qu. Median Mean 3rd Qu. Max. 1.00 3.00 4.00 13.86 24.50 47.00 > summary(C) Min. 1st Qu. Median Mean 3rd Qu. Max. 1.00 5.00 11.00 11.09 16.00 20.00 グラフ化なしで統計処理をするな、という例だな。
67 :卵の名無しさん :2020/03/17(火) 00:15:40.58 ID:+SeXrvpf.net 中央値比較での統計悪用シミュレーションプログラム sim <- function(){ A=c(sort(sample(1:4,17,rep=T)),4,sort(sample(5:50,17,rep=T))) C=c(sort(sample(1:11,23,rep=T)),11,sort(sample(11:20,23,rep=T))) print(summary(A)) print(summary(C)) par(mfrow=c(2,1)) plot(table(A),bty='l',xlim=c(1,50),col='red') plot(table(C),bty='l',xlim=c(1,50),col='blue') } sim()
68 :卵の名無しさん :2020/03/17(火) 00:19:03.10 ID:+SeXrvpf.net こんな感じ https://i.imgur.com/kL3pnn3.png https://i.imgur.com/3HYkd5O.png
69 :卵の名無しさん :2020/03/17(火) 11:29:17.69 ID:qm/RkVMU.net >>64 ド底辺頭脳はサイコロを6回投げたら1の目が1回でているとう風な計算しかできんみたいだな。 信頼区間とか確率分布という概念での思考ができなんよね。 周辺度数固定して二項分布で乱数発生させてχ二乗検定したときのp値の分布は以下のようになる。 https://i.imgur.com/5djScD1.png https://i.imgur.com/RsFr7o4.png 危険率1%で有効といえる確率はほぼ5割であることがわかる。
70 :卵の名無しさん :2020/03/17(火) 11:38:17.66 ID:qm/RkVMU.net >>60 エントリーに5以下の数値があるとχ二乗検定だとYatesの補正を使うけど、この程度の標本数は臨床では普通だよ。 中国のアビガンの予備試験も標本数が実薬群ととコントロール群を合わせて80だぞ。 足りないのはオマエの学力だよ。
71 :卵の名無しさん :2020/03/17(火) 12:02:20.58 ID:RqUjOCLx.net >>70 まあ 学園カースト最下位のFランさんが必死ですわね 昨日は先輩方からも国試合格のご報告がたくさんありましたが Fランさんのお知り合いの方たちはいかがですかぁ?W by 都内Sラン女子高生
72 :卵の名無しさん :2020/03/17(火) 12:41:00.33 ID:+SeXrvpf.net >>71 答えられない低学力を隠そうと必死だな。 どうやって解くかわからんの?
73 :卵の名無しさん :2020/03/17(火) 19:55:56.01 ID:RqUjOCLx.net >>72 わかりました 犯人はアナタですね Fランさん
74 :卵の名無しさん :2020/03/17(火) 21:23:49.50 ID:qm/RkVMU.net ようやく完成した # 球面上の一様分布 vertex <- function(r=1){ x=runif(1,-1,1) # x ~ 一様分布[-1,1] phi=runif(1,-pi,pi) # φ ~ 一様分布[-π,π] y=sqrt(1-x^2)*cos(phi) # √(1-x^2)*cos(φ) z=sqrt(1-x^2)*sin(phi) # √(1-x^2)*sin(φ) r*c(x,y,z) } https://i.imgur.com/xX0mTim.png https://i.imgur.com/H7hs9w8.png # 凸包四面体の体積 sim <- function(r=1,print=F){ v4=replicate(4,vertex()) # 4点の直交座標 if(print) print(v4) abs(det(v4[,2:4]-v4[,1]))/6 # 四面体の体積 } sim(print=T) k=1e5 tetra=replicate(k,sim()) # k回のシミュレーション mean(tetra) summary(tetra) BEST::plotPost(tetra)
75 :卵の名無しさん :2020/03/17(火) 21:46:10.60 ID:qm/RkVMU.net # 単位球体内部一様分布の4点で四面体をつくるとき rm(list=ls()) vertex <- function(){ x=runif(1,-1,1) # x ~ 一様分布[-1,1] phi=runif(1,-pi,pi) # φ ~ 一様分布[-π,π] r=runif(1,0,1) # r ~ 一様分布[0,1] x=r^(1/3)*x # r^(1/3)*x y=r^(1/3)*sqrt(1-x^2)*cos(phi) # r^(1/3)√(1-x^2)*cos(φ) z=r^(1/3)*sqrt(1-x^2)*sin(phi) # r^(1/3)√(1-x^2)*sin(φ) c(x,y,z) } vtx=replicate(5000,vertex()) par(mfrow=c(3,1)) x=vtx[1,] ; hist(x,col='pink') y=vtx[2,] ; hist(y,col='orange') z=vtx[3,] ; hist(z,col='darkgreen') rgl::plot3d(x,y,z, col="skyblue") par(mfrow=c(1,1)) # 四面体の体積 sim <- function(r=1,print=F){ v4=replicate(4,vertex()) # 4点の直交座標 if(print) print(v4) abs(det(v4[,2:4]-v4[,1]))/6 # 四面体の体積 } sim(print=T) k=1e5 tetra=replicate(k,sim()) # k回のシミュレーション mean(tetra) summary(tetra) BEST::plotPost(tetra)
76 :卵の名無しさん :2020/03/17(火) 21:47:20.28 ID:qm/RkVMU.net core i7 Memory 16G はシミュレーションが捗って( ・∀・)イイ!!
77 :卵の名無しさん :2020/03/20(金) 07:04:24.48 ID:LvYpiZ17.net 臨床医学は確率事象を扱う稼業なので統計処理は必須。 これができない医師はアホといえる。 某国の新型コロナ感染症の有病率を0.1、 PCR検査の感度を0.7 特異度を0.9とする 検査陽性陰性の人を無作為に100人ずつ集めて、カレーを頻回に食べているかを調査した結果が 以下の通りであった。 検査陽性 検査陰性 カレー頻食 a(=36) b(=60) カレー稀食 c(=64) d(=40) カレーを頻食すると新型コロナ感染に罹りにくいと結論できるか?
78 :卵の名無しさん :2020/03/21(土) 10:28:11.17 ID:eZCI10Rh.net Parameter Total (n=83) Severe/critica l Group (n=25) Ordinary Group (n=58) P Value Number of lobes 5 (4-5) 5 (5-5) 5 (3-5) 0.003
79 :卵の名無しさん :2020/03/21(土) 10:29:01.74 ID:eZCI10Rh.net median (interquartile range, IQR)
80 :卵の名無しさん :2020/03/21(土) 17:41:24.46 ID:o3rnHhhm.net library(coin) # coin::wilcox_test(c(Ea,La) ~ factor(c(rep('Ea',n1),rep('La',n2)))) f4 <- function(x){ # filer p.value=0.003 # pv=t.test(sc[x[1],],or[x[2],])$p.value SC=sc[x[1],] OR=or[x[2],] w=coin::wilcox_test(c(SC,OR) ~ factor(c(rep('SC',25),rep('OR',58))),distribution="exact") pv=coin::pvalue(w) if(round(pv,3)==0.003) return(c(x[1],x[2])) else(return(c(0,0))) }
81 :卵の名無しさん :2020/03/22(日) 09:18:45.02 ID:/FOXgWJQ.net PCR <- function(n=1000,r=10,sen=0.7,spc=0.9,k=1e6,print=TRUE){ # n:検査件数 r:検査陽性数 sen=感度 spc=特異度 k:発生乱数の数,print:グラフ描画 prev=rbeta(k,1+r,1+n-r) # 有病率の事前分布を一様分布と仮定 PPV=prev*sen/( prev*sen + (1-prev)*(1-spc) ) # 検査特性を加味してPPV計算 adjtd.prev = PPV*(1-spc)/(sen - PPV*(sen+spc-1)) # そのPPVから有病率を計算 # NPVでも同じ結果 # NPV=(1-prev)*spc/( (1-prev)*spc + prev*(1-sen)) # 検査特性を加味してNPV計算 # adjtd.prev = (1-NPV)*spc/(spc - NPV*(sen+spc-1)) # そのNPVから有病率を計算 mean=mean(adjtd.prev) # 期待値 median=median(adjtd.prev) # 中央値 mode=density(adjtd.prev)$x[which.max(density(adjtd.prev)$y)] # 最頻値 LU=unlist(HDInterval::hdi(adjtd.prev))[1:2] # 95%信頼区間 re=c(mean=mean,median=median,mode=mode,LU) # 事後有病率の代表値 if(print){ # ヒストグラム描画 par(mfrow=c(2,2)) hist(prev,freq=F,xlim=c(0,0.5), main='prevalence') ; lines(density(prev)) hist(PPV,freq=F,xlim=c(0,0.5),main='PPV') ; lines(density(PPV)) hist(adjtd.prev,freq=F,xlim=c(0,0.5),main='adjusted prevalence ') ; lines(density(adjtd.prev)) BEST::plotPost(adjtd.prev) ; lines(density(adjtd.prev),col='skyblue') par(mfrow=c(1,1)) } print(round(re,4)) # 概算値表示 invisible(list(re,adjtd.prev)) # 代表値と事後有病率を非表示で返す } PCR(100,1) PCR(1000,10) PCR(10000,100) PCR(100000,1000)
82 :卵の名無しさん :2020/03/22(日) 20:13:20.91 ID:/FOXgWJQ.net # stanでMCMCさせることにした。 data{ // corona.stan int n; // 1000 int x; // 10 real<lower=0,upper=1> sen; // 0.7 real<lower=0,upper=1> spc; // 0.9 } parameters{ real<lower=0,upper=1> prev; // prevalence } transformed parameters{ real<lower=0,upper=1> p; p = prev*sen + (1-prev)*(1-spc) ; // probability of positive test result } model{ x ~ binomial(n,p); prev ~ beta(1,1); // prev ~ uniform(0,1) }
83 :卵の名無しさん :2020/03/22(日) 20:14:04.17 ID:/FOXgWJQ.net stanはコンパイルに時間がかかって面倒。 jagsの方が早くて( ・∀・)イイ!! rm(list=ls()) graphics.off() options(scipen = 4) options(digits=5) library(rjags) N=1000 ; X=10 dataList=list(n=N,x=X,sen=0.7,spc=0.9) modelstring=paste0(' model { prev ~ dbeta(1,1) p <- prev*sen + (1-prev)*(1-spc) x ~ dbin(p,n) } ') writeLines(modelstring,'TEMPmodel.txt') jagsModel = jags.model( file="TEMPmodel.txt" ,data=dataList) update(jagsModel) codaSamples = coda.samples( jagsModel , variable=c("prev","p"), n.iter=50000 ) js=as.matrix(codaSamples) BEST::plotPost(js[,'prev']) round(HDInterval::hdi(js[,'prev'])[1:2],7)
84 :卵の名無しさん :2020/03/23(月) 12:03:19.86 ID:a9wYdNxH.net 何の責任もとれないド底辺ほど 他人の責任にうるさく自分の権利ばかりしつこく要求しがちですわね コンプ薬屋のことですけど Fランのアナタも同類ですわよね 早く呼んでくればぁ by 都内Sラン女子高生
85 :卵の名無しさん :2020/03/23(月) 13:38:33.35 ID:a9wYdNxH.net >>1 よそのスレで騒いでるみたいですけど Fランのアナタには知性があるとは 医師の方々は考えてはいないようですよ by 都内Sラン女子高生
86 :卵の名無しさん :2020/03/23(月) 14:45:03.80 ID:a9wYdNxH.net >>1 国試に行き詰まって統計遊びですか こういうときハゲワロスって言うんでしたっけ? by ドSなSラン女子高生
87 :卵の名無しさん :2020/03/23(月) 14:48:48.41 ID:a9wYdNxH.net タバコを吸わないものは自室に禁煙の張り紙など貼らない訳ですが 裏口裏口と盛んにワメキ回るアナタには 身に覚えがあると言うことでよろしいですねW 容疑者はオマエだぁっ(コナ○口調で) by ドSなSラン女子高生
88 :卵の名無しさん :2020/03/23(月) 15:00:20.99 ID:a9wYdNxH.net アナタの期待値は100発0中だぁっ ち○ぴょろすぽ○ん by ドSなSラン女子高生
89 :卵の名無しさん :2020/03/23(月) 18:19:32.21 ID:a9wYdNxH.net 自分で作った統計ごっこのスレがここにあるのに なーんでよそのスレに迷惑をかけに行くのかなん やっぱりFラン脳のゆえんなのかしらねんのねん by ドSなSラン女子高生
90 :卵の名無しさん :2020/03/25(水) 09:16:02.76 ID:htH1Knlx.net >>88 これやってみ、ド底辺頭脳のオマエには無理だろうな。 24時間待ってやるぞ。 ドツボ13は100発0中 ドツボ14は10発0中 ドツボ15は1発0中 とする。 各々10000発撃ったときドツボの命中数の期待値はいくらか?
91 :卵の名無しさん :2020/03/25(水) 09:16:47.51 ID:htH1Knlx.net 正確な感度も特異度もわからないから有病率は計算できないという主張が数学板でされたから、 感度・特異度を確率変数としてベイズ階層モデルを組めば計算できること投稿しておいた。 最頻値と標準偏差からベータ分布のパラメータを計算するのに連立方程式を解く方が面倒だった。
92 :卵の名無しさん :2020/03/25(水) 14:59:53.46 ID:Jsfh+Odg.net きゃー、統計ごっこのクルクルパー来たーっW (キモワロス) ド底辺仲間のコンプ薬屋はまだですか? by ドSなSラン女子高生
93 :卵の名無しさん :2020/03/25(水) 15:55:21.10 ID:C3cAic+l.net >>92 >>88 これやってみ、ド底辺頭脳のオマエには無理だろうな。 24時間待ってやるぞ。 ドツボ13は100発0中 ドツボ14は10発0中 ドツボ15は1発0中 とする。 各々10000発撃ったときドツボの命中数の期待値はいくらか?
94 :卵の名無しさん :2020/03/26(木) 10:22:15.80 ID:F34mrbtj.net >>93 国試に合格できなかった還暦ジジイって本当ですかぁ? アチコチで反日レスしてるみたいですねー by ドSなSラン女子高生
95 :卵の名無しさん :2020/03/26(木) 10:47:49.20 ID:Z3UzEFJB.net >>94 ドツボの宿題まだぁ?
96 :卵の名無しさん :2020/03/26(木) 11:18:27.82 ID:F34mrbtj.net >>95 Fラン高の宿題はFランのアナタ自身でされてくださいね Sランのアタシはアナタと違ってヒマではございませんのでW by ドSなSラン女子高生
97 :卵の名無しさん :2020/03/26(木) 19:27:21.07 ID:O1C2WieZ.net >>96 ドツボ問題もできないドアホ。 これやってみ、ド底辺頭脳のオマエには無理だろうな。 24時間待ってやるぞ。 ドツボ13は100発0中 ドツボ14は10発0中 ドツボ15は1発0中 とする。 各々10000発撃ったときドツボの命中数の期待値はいくらか?
98 :卵の名無しさん :2020/03/27(金) 07:08:38.46 ID:BC/LnmbB.net 95人調べて47人陽性。感度50-70% 特異度90%前後として この集団の有病率はstanでMCMCして計算すると mean lower upper 0.7558614 0.5214634 0.9999897 感度60% 特異度90%に固定すると > pn2ip(95,47) infected prevalence 75.0000000 0.7894737
99 :卵の名無しさん :2020/04/04(土) 05:13:53.27 ID:ApKEo+4Y.net 麻薬取締法違反 札幌の薬剤師と法人を略式起訴 病院で管理する医療用麻薬の数量について道に虚偽の届け出をしたとして、 札幌区検は19日までに麻薬取締法違反の罪で、札幌ひばりが丘病院(札幌 市厚別区)に勤務していた30代の薬剤師の男と、同病院を運営する医療法人 潤和会(同区)を略式起訴した。札幌簡裁がそれぞれ20万円以下の罰金刑 を言い渡す見通し。略式命令請求によると、男は2015年11月、道に対し、 院内で使用した麻薬の使用量を偽って届け出るなどしたとされる。潤和会には 法人も罰する同法の両罰規定を適用した。 出典:北海道新聞 平成30年10月19日付
100 :卵の名無しさん :2020/04/07(火) 08:43:20.38 ID:je8pKp8V.net 晋型コロナ肺炎に感度0.9,特異度0.9の迅速検査が開発されたと仮定する。 日本人1億2595万人からX人を無作為抽出して有病率を推定したい。 有病率の99%信頼区間幅を1%以内で検定したい。 何人を抽出すれば十分といえるか? ss2sbj <- function(SEN,SPC,range=0.01,CONF=0.99,print=FALSE){ pr2sbj <- function(p0,sen=SEN,spc=SPC,R=range,conf.level=CONF){ z = qnorm(1-(1-conf.level)/2) p = p0*sen+(1-p0)*spc n = 4*z^2*p*(1-p)/R^2 return(n) } if(print){ prevalence=seq(0,1,by=0.001) plot(prevalence,sapply(prevalence,pr2sbj), ylab='subjects',bty='l',type='l',lwd=2)} optimize(function(x) pr2sbj(x,SEN,SPC), c(0,1) ,maximum = TRUE) }
101 :卵の名無しさん :2020/04/22(水) 06:24:37.56 ID:XcybEEeM.net "n(=100)人の中から無作為にm(=10)人選んだらその中に少なくとも一人の感染者がいた。 全体で何人の感染者がいるかの期待値と95%信頼区間を求めよ。" rm(list=ls()) fn <- function(n,m,print=FALSE){ library(gmp) # 離散量で計算 x=0:n # 感染者数:x, 非感染数:n-x pmf=1- chooseZ(n-x,m)/chooseZ(n,m) # 1 - (m人全員非感染の確率) pdf=pmf/sum(pmf) # 確率密度関数化して Ev=sum(x*pdf) # 期待値を計算 if(print) print(Ev) E=as.numeric(Ev) # 信頼区間計算計算のために連続量とする y=seq(0,n,by=0.05) # infected p=1- choose(n-y,m)/choose(n,m) # Pr[at least one infected among m] d=p/sum(p) # pmf -> pdf (e=sum(y*d)) # mean of probability if(print){plot(y,d,bty='l',xlab='infected',ylab='density',type='l') print(quantile(n*p,prob=c(0.025,0.05,0.5,0.95,0.975))) print(HDInterval::hdi(n*p)[1:2])} return(E) } fn(100,5,print=T) n=100 m=1:n plot(m,sapply(m,function(x) fn(n,x)), type='p',pch=19,bty='l',ylim=c(0,100), xlab='有名人',ylab='期待値',main='総人口:100人')
145 KB
新着レスの表示
掲示板に戻る
全部
前100
次100
最新50
read.cgi ver.24052200
本文 スレッドタイトル 投稿者