臨床統計もおもしろいですよ、その3

1 ：卵の名無しさん：2020/03/05(木) 20:17:05.03 ID:naSB8128.net

　
　内科認定医受験の最低限の知識、
　製薬会社の示してくる臨床データ、
　論文の考察、
　論文を書くときの正当性、
　というのが、臨床統計の今までの目的の大きい部分でしたが、
　
　AI=機械学習の基本も、結局は統計学と確率に支配されます。
　そういう雑多な話をするスレです。
　
※前スレ
臨床統計もおもしろいですよ、その１
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人の集団（クラスター）で感染者は１人として１日に延べ１０回・人と接触し、
１人１回あたりの感染確率を１％、感染期間30日、潜伏期5日として
SEIRモデルで計算すると、感染のピークは110日めでとても１〜２週間が山場とは言えない。


> 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
受験を控えている高校生です

他のスレでアナタはコンプ性Ｆラン症候群という不治の病であるとお聞きしました
高学歴の私と違ってＦランであるアナタは今はどのようなお仕事をされてるんですか？

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モデルで有病率を１％に固定して、集団のサイズを変化させてシミュレーションしてみたけどピークは変わらないな。
このモデルでは集会規模の大小には影響されないということになるな。

https://i.imgur.com/343K91V.png

有病率を変化させて流行の変遷をグラフにすると、

https://i.imgur.com/SZ15LKT.png
https://i.imgur.com/gnJVFnd.png

有病率を４０％くらいに引き上げるとオリンピックのときには流行が収束していることになるｗｗｗｗ

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
同じアルゴリズムを使って有病率が１/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

"
東京都は２１日、新型コロナウイルスの感染拡大を防ぐため、２２日から３月１５日までの３週間、都が主催する５００人以上の大規模な屋内イベントは、原則延期か中止にする方針を発表した.
"
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) # ｎの順列
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 # １の目の出た回数
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してｐ値の信頼区間を出してみた。
#
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

ｐ値の期待値は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

ｆｎ　＜−　ｆｕｎｃｔｉｏｎ（ｎ）｛
ｌｉｂｒａｒｙ（ｇｍｐ）
ｇｒ＝ｅｖａｌ（ｓｔｒ２ｌａｎｇ（ｐａｓｔｅ（�Wｅｘｐａｎｄ．ｇｒｉｄ（�W，ｐａｓｔｅ（ｒｅｐ（�W１：６�W，ｎ），ｃｏｌｌａｐｓｅ＝�V，�V），�W）�W）））
ｒ＝ａｓ．ｂｉｇｑ（ｓｕｍ（ａｐｐｌｙ（ｇｒ，１，ｆｕｎｃｔｉｏｎ（ｘ）　ｐｒｏｄ．ｂｉｇｚ（ｘ）％％６＝＝０）））／ａｓ．ｂｉｇｑ（ｎｒｏｗ（ｇｒ））
ｒ２＝ｃａｐｔｕｒｅ．ｏｕｔｐｕｔ（ｒ）［２］
ｓｕｂｓｔｒ（ｒ２，５，ｎｃｈａｒ（ｒ２））
｝
ｆｏｒ（ｉ　ｉｎ　１：９）　ｃａｔ（ｉ，�V：�V，ｆｎ（ｉ），�V\ｎ�V）

＃　ｓｔｒ２ｌａｎｇ　ｓｔｒ２ｅｘｐｒｅｓｓｉｏｎ
？ｓｔｒ２ｌａｎｇ
ｒｅｐ（�W１：６�W，５）
ｐａｓｔｅ（ｒｅｐ（�W１：６�W，５），ｃｏｌｌａｐｓｅ＝�V，�V）
ｓｔｒ２ｌａｎｇ（ｐａｓｔｅ（�Wｅｘｐａｎｄ．ｇｒｉｄ（�W，ｐａｓｔｅ（ｒｅｐ（�W１：６�W，５），ｃｏｌｌａｐｓｅ＝�V，�V），�V）�V））
ｅｖａｌ（ｓｔｒ２ｌａｎｇ（ｐａｓｔｅ（�Wｅｘｐａｎｄ．ｇｒｉｄ（�W，ｐａｓｔｅ（ｒｅｐ（�W１：６�W，５），ｃｏｌｌａｐｓｅ＝�V，�V），�V）�V）））

（ｓ＝ｒｅｐ（�W１：６�W，５））
（ｓｔｒ＝ｐａｓｔｅ（ｓ，ｃｏｌｌａｐｓｅ＝�V，�V））
（ｌａｎｇ＝ｓｔｒ２ｌａｎｇ（ｐａｓｔｅ（�Wｅｘｐａｎｄ．ｇｒｉｄ（�W，ｓｔｒ，�V）�V）））
ｅｖａｌ（ｌａｎｇ）

41 ：卵の名無しさん：2020/03/15(日) 16:22:13.79 ID:5ItQZoKJ.net

>>1
受験を控えている女子高校生です

他のスレでアナタはコンプ性Ｆラン症候群という不治の病であるとお聞きしました
高学歴の私と違ってＦランであるアナタは今はどのようなお仕事をされてるんですか？

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 都内高学歴ＪＫ

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
早くＦラン仲間のコンプ薬屋を呼んで来なよ
お勉強してくれるって言ってるお

by 都内高学歴ＪＫ

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 日あたり新たに新型コロナウイルス感染症を疑って外来を受診する患者数）＝
　　　（0-14 歳人口）×0.18／100＋（15-64 歳人口） ×0.29／100＋（65 歳以上人口） ×0.51／100

（２）（ピーク時において1 日あたり新型コロナウイルス感染症で入院治療が必要な患者数）＝
　　　（0-14 歳人口）×0.05／100＋（15-64 歳人口）×0.02／ 100＋（65 歳以上人口） ×0.56／100

（３）（ピーク時において1 日あたり新型コロナウイルス感染症で重症者として治療が必要な患者数）＝
　　　（0-14 歳人口）×0.002／100＋（15-64 歳人口） ×0.001／100＋（65 歳以上人口） ×0.018／100

注１）ピーク時は、各都道府県等において疫学的関連性が把握できない程度に感染が拡大した時点から概ね３か月後に到来すると推計されている。ただし、公衆衛生上の対策を行うことにより、ピークが下がるとともに後ろ倒しされる。
注２）重症者とは、集中治療や人工呼吸器を要する管理が必要な患者を指す。
注３）当該計算式は、都道府県等の単位以下における医療提供体制を確保するためのものであるとともに、各都道府県等によってピークを迎える時期が異なるため、
全国の人口を用いて計算することや単純に各自治体が算出するピークの数値を足し合わせることは、不適切な取扱いとなることに留意いただきたい。なお、当該計算式については、
今後新たな知見等により変更される可能性がある。
注４）実際には、ピーク時に至るまでの日々の患者数の増加はばらつきがあり、増加曲線は推計通りの形にならない可能性が高いため、
現実の患者の発生動向も踏まえて適切に体制を確保することが必要。
注５）当該計算式については、今後新たな知見等により変更される可能性がある。

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