臨床統計もおもしろいですよ、その3
内科認定医受験の最低限の知識、
製薬会社の示してくる臨床データ、
論文の考察、
論文を書くときの正当性、
というのが、臨床統計の今までの目的の大きい部分でしたが、
AI=機械学習の基本も、結局は統計学と確率に支配されます。
そういう雑多な話をするスレです。
※前スレ
臨床統計もおもしろいですよ、その１
https://egg.2ch.sc/test/read.cgi/hosp/1493809494/
臨床統計もおもしろいですよ、その2
https://egg.2ch.sc/test/read.cgi/hosp/1540905566/
１人１回あたりの感染確率を１％、感染期間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
https://egg.2ch.sc/test/read.cgi/hosp/1529634250/
古根高（サッポロファクトリーメンタルクリニック院長）（五浪で底辺医大に進学した知恵遅れ）のご尊顔
https://i.imgur.com/u8rBskz.jpg
トラブルメーカーゆえ市立札幌病院を解雇された過去あり
https://i.imgur.com/upWjhiA.jpg
nCov2019 for studying COVID-19 coronavirus outbreak
https://guangchuangyu.github.io/nCov2019/
まず、Rを最新版に更新して
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("BiocStyle")
remotes::install_github("GuangchuangYu/nCov2019", dependencies = TRUE)
が必要だった。
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
受験を控えている高校生です
他のスレでアナタはコンプ性Ｆラン症候群という不治の病であるとお聞きしました
高学歴の私と違ってＦランであるアナタは今はどのようなお仕事をされてるんですか？
http://www.bcloud.org/e/
これは使えるな
http://www.bcloud.org/e/session/a9b2a4a3e4ef04951a3e45743536a8f3/download/dataDownloadWorld?w=
このモデルでは集会規模の大小には影響されないということになるな。
https://i.imgur.com/343K91V.png
有病率を変化させて流行の変遷をグラフにすると、
https://i.imgur.com/SZ15LKT.png
https://i.imgur.com/gnJVFnd.png
有病率を４０％くらいに引き上げるとオリンピックのときには流行が収束していることになるｗｗｗｗ
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)
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)
東京都は２１日、新型コロナウイルスの感染拡大を防ぐため、２２日から３月１５日までの３週間、都が主催する５００人以上の大規模な屋内イベントは、原則延期か中止にする方針を発表した.
"
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='感染拡大確率')
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
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)
}
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
これバグがあるな
バグ修正版
https://rio2016.2ch.sc/test/read.cgi/math/1501755792/789
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))
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)
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)
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)
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)
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
}
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')
何回目に賭けるのがベストか？"
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')
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)
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);
}
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
誰が正規入学かを確かめる必要がある。審議官の女医（嘘つきかどうかは不明）が無作為に新入生を呼び出して
「あなたのいうことが正しければ手コキかRをしてあげる、そうでなければセンズリを命じる」という試問をする。
論理的思考ができる正規入学生はこの質問に回答してRをしてもらっているが、裏口入学生は正解がだせずにセンズリを命じられている。
正規入学生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
検査陽性 検査陰性
カレー頻食 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()
検査陽性陰性の人を無作為に50人ずつ集めて、カレーを頻回に食べているかを調査した結果が
以下の通りであった。
検査陽性 検査陰性
カレー頻食 a(=18) b(=30)
カレー稀食 c(=32) d(=20)
カレーを頻食すると新型コロナ感染に罹りにくいと結論できるか？ 危険率0.05で検定せよ。
有病率によって有意かどうか判断が変わるのは興味ぶかい。
https://i.imgur.com/qjjIjnm.png
周辺度数を固定して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割という結果がえられた。
タイプミス修正
症数例から
↓
少数例
で知った
例
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)
ｌｉｂｒａｒｙ（ｇｍｐ）
ｇｒ＝ｅｖａｌ（ｓｔｒ２ｌａｎｇ（ｐａｓｔｅ（?ｅｘｐａｎｄ．ｇｒｉｄ（?，ｐａｓｔｅ（ｒｅｐ（?１：６?，ｎ），ｃｏｌｌａｐｓｅ＝?，?），?）?）））
ｒ＝ａｓ．ｂｉｇｑ（ｓｕｍ（ａｐｐｌｙ（ｇｒ，１，ｆｕｎｃｔｉｏｎ（ｘ） ｐｒｏｄ．ｂｉｇｚ（ｘ）％％６＝＝０）））／ａｓ．ｂｉｇｑ（ｎｒｏｗ（ｇｒ））
ｒ２＝ｃａｐｔｕｒｅ．ｏｕｔｐｕｔ（ｒ）［２］
ｓｕｂｓｔｒ（ｒ２，５，ｎｃｈａｒ（ｒ２））
｝
ｆｏｒ（ｉ ｉｎ １：９） ｃａｔ（ｉ，?：?，ｆｎ（ｉ），?\ｎ?）
＃ ｓｔｒ２ｌａｎｇ ｓｔｒ２ｅｘｐｒｅｓｓｉｏｎ
？ｓｔｒ２ｌａｎｇ
ｒｅｐ（?１：６?，５）
ｐａｓｔｅ（ｒｅｐ（?１：６?，５），ｃｏｌｌａｐｓｅ＝?，?）
ｓｔｒ２ｌａｎｇ（ｐａｓｔｅ（?ｅｘｐａｎｄ．ｇｒｉｄ（?，ｐａｓｔｅ（ｒｅｐ（?１：６?，５），ｃｏｌｌａｐｓｅ＝?，?），?）?））
ｅｖａｌ（ｓｔｒ２ｌａｎｇ（ｐａｓｔｅ（?ｅｘｐａｎｄ．ｇｒｉｄ（?，ｐａｓｔｅ（ｒｅｐ（?１：６?，５），ｃｏｌｌａｐｓｅ＝?，?），?）?）））
（ｓ＝ｒｅｐ（?１：６?，５））
（ｓｔｒ＝ｐａｓｔｅ（ｓ，ｃｏｌｌａｐｓｅ＝?，?））
（ｌａｎｇ＝ｓｔｒ２ｌａｎｇ（ｐａｓｔｅ（?ｅｘｐａｎｄ．ｇｒｉｄ（?，ｓｔｒ，?）?）））
ｅｖａｌ（ｌａｎｇ）
受験を控えている女子高校生です
他のスレでアナタはコンプ性Ｆラン症候群という不治の病であるとお聞きしました
高学歴の私と違ってＦランであるアナタは今はどのようなお仕事をされてるんですか？
高学歴ならこれに答えてみ！
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)
# 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)
# 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)
ド底辺頭脳でも数くらい数えられるんだろ？
1つのサイコロを10000回投げ、出たの目の積が6の倍数になる確率を分数でもとめよ。
どうでもいいから早くコンプ薬屋を呼んで来なよ
そんなだからモテないんだお
by 都内高学歴ＪＫ
数も数えられんの？ガイジ？
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)
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)
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)
早くＦラン仲間のコンプ薬屋を呼んで来なよ
お勉強してくれるって言ってるお
by 都内高学歴ＪＫ
ド底辺頭脳でも数くらい数えられるんだろ？
1つのサイコロを10000回投げ、出たの目の積が6の倍数になる確率を分数でもとめよ。
今後、国内で新型コロナウイルス感染症患者数が大幅に増えたときに備え、各都道府県、保健所設置市及び特別区（以下「都道府県等」という。）における外来を受診する患者数等について、
以下の数式を用いて計算いただき、ピーク時の医療需要の目安としてご活用の上、必要な医療提供体制を確保していただくようお願いいたします。
（１）（ピーク時において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
注１）ピーク時は、各都道府県等において疫学的関連性が把握できない程度に感染が拡大した時点から概ね３か月後に到来すると推計されている。ただし、公衆衛生上の対策を行うことにより、ピークが下がるとともに後ろ倒しされる。
注２）重症者とは、集中治療や人工呼吸器を要する管理が必要な患者を指す。
注３）当該計算式は、都道府県等の単位以下における医療提供体制を確保するためのものであるとともに、各都道府県等によってピークを迎える時期が異なるため、
全国の人口を用いて計算することや単純に各自治体が算出するピークの数値を足し合わせることは、不適切な取扱いとなることに留意いただきたい。なお、当該計算式については、
今後新たな知見等により変更される可能性がある。
注４）実際には、ピーク時に至るまでの日々の患者数の増加はばらつきがあり、増加曲線は推計通りの形にならない可能性が高いため、
現実の患者の発生動向も踏まえて適切に体制を確保することが必要。
注５）当該計算式については、今後新たな知見等により変更される可能性がある。
0.05,0.02,0.56,
0.002,0.001,0.018)
(mat=matrix(v*100,3,byrow=T))
童貞仲間でしょぉお
by 都内高学歴ＪＫ
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人の死者予想。
コンプはこれに答がだせなかったんだが、あんたは出せる？
解答できなきゃ、コンプと同レベルの頭脳と認定してあげよう。
24時間待ってあげる。
問題
某国の新型コロナ感染症の有病率を0.1、PCR検査の感度を0.7 特異度を0.9とする
検査陽性陰性の人を無作為に50人ずつ集めて、カレーを頻回に食べているかを調査した結果が
以下の通りであった。
検査陽性 検査陰性
カレー頻食 a(=18) b(=30)
カレー稀食 c(=32) d(=20)
カレーを頻食すると新型コロナ感染に罹りにくいと結論できるか？ 危険率0.05で検定せよ。
今後、国内で新型コロナウイルス感染症患者数が大幅に増えたときに備え、各都道府県、保健所設置市及び特別区（以下「都道府県等」という。）における外来を受診する患者数等について、
以下の数式を用いて計算いただき、ピーク時の医療需要の目安としてご活用の上、必要な医療提供体制を確保していただくようお願いいたします。
（１）（ピーク時において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
問題 ： ある都市でピーク時に外来1000人、入院600人、重症20人であったとすると、この都市の14歳以下の人口は何人と推測されるか？
評価するには症例数が２桁少ないですわね
それはさておき昔の教授選では積み上げたペーパーの高さを比べたと聞きますが
Ｆランのアナタは採択されたペーパーの１本もあるのかしら？Ｗ
by 都内Ｓラン女子高生
足りないのは標本数でなくてオマエの頭だよ。
国家緊急予防および制御薬物工学技術研究センターと深セン第三人民病院LeレイおよびLi英Yチームによって完成された「ファピラビル」。
「新しいコロナウイルス肺炎（COVID-19）患者の安全性と有効性に関する臨床研究」は、予備的な臨床結果に達しています（登録番号：
ChiCTR2000029600）。
研究により、ファビピラビルは、ウイルスクリアランスを促進することにより、新しいコロナウイルス肺炎の進行を緩和することが示されています。研究結果は、中国工学院のジ
ャーナルに提出されました
合計80人の患者がこの臨床試験に登録され、35人の患者がファピラビルで治療され、対照群は年齢、性別、および疾患の重症度が治療群と一致した新冠動脈肺炎の患者45人
でした。 iタブレット治療。薬物投与からウイルス除去までの時間の中央値、治療の14日目の胸部画像の改善率、および安全性を2つのグループ間で比較しました。
結果は、ウイルス除去の観点から、治療後のファピラビル試験群の患者のウイルス除去時間の中央値（ウイルス核酸陰性）は、それぞれ4日と11日であった対照群のそれよりも
有意に短いことを示しました。別の重要な指標として、患者の胸部画像の改善、治療群と対照群の改善率はそれぞれ91.43％と62.22％でした。同時に、対照群と比較し
て、Rビル治療群は副作用が少なく、忍容性が良好でした。 2つの患者グループのベースライン特性はすべて同等でした。
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
https://i.imgur.com/3uZNIwI.png
PCR陰転までの日数の中央値が
アビガン群３５例で４日、対照群で１１日としても
次のような分布なら有効とは言い難い。
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
グラフ化なしで統計処理をするな、という例だな。
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()
https://i.imgur.com/kL3pnn3.png
https://i.imgur.com/3HYkd5O.png
ド底辺頭脳はサイコロを６回投げたら１の目が１回でているとう風な計算しかできんみたいだな。
信頼区間とか確率分布という概念での思考ができなんよね。
周辺度数固定して二項分布で乱数発生させてχ二乗検定したときのｐ値の分布は以下のようになる。
https://i.imgur.com/5djScD1.png
https://i.imgur.com/RsFr7o4.png
危険率１％で有効といえる確率はほぼ５割であることがわかる。
エントリーに５以下の数値があるとχ二乗検定だとYatesの補正を使うけど、この程度の標本数は臨床では普通だよ。
中国のアビガンの予備試験も標本数が実薬群ととコントロール群を合わせて80だぞ。
足りないのはオマエの学力だよ。
まあ 学園カースト最下位のＦランさんが必死ですわね
昨日は先輩方からも国試合格のご報告がたくさんありましたが
Ｆランさんのお知り合いの方たちはいかがですかぁ？Ｗ
by 都内Ｓラン女子高生
答えられない低学力を隠そうと必死だな。
どうやって解くかわからんの？
わかりました
犯人はアナタですね Ｆランさん
# 球面上の一様分布
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()) # ４点の直交座標
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)
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()) # ４点の直交座標
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)
これができない医師はアホといえる。
某国の新型コロナ感染症の有病率を0.1、
PCR検査の感度を0.7 特異度を0.9とする
検査陽性陰性の人を無作為に100人ずつ集めて、カレーを頻回に食べているかを調査した結果が
以下の通りであった。
検査陽性 検査陰性
カレー頻食 a(=36) b(=60)
カレー稀食 c(=64) d(=40)
カレーを頻食すると新型コロナ感染に罹りにくいと結論できるか？
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
# 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)))
}
# 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)
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)
}
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)
他人の責任にうるさく自分の権利ばかりしつこく要求しがちですわね
コンプ薬屋のことですけど
Fランのアナタも同類ですわよね
早く呼んでくればぁ
by 都内Ｓラン女子高生
よそのスレで騒いでるみたいですけど
Fランのアナタには知性があるとは
医師の方々は考えてはいないようですよ
by 都内Ｓラン女子高生
国試に行き詰まって統計遊びですか
こういうときハゲワロスって言うんでしたっけ？
by ドＳなＳラン女子高生
裏口裏口と盛んにワメキ回るアナタには
身に覚えがあると言うことでよろしいですねＷ
容疑者はオマエだぁっ（コナ○口調で）
by ドＳなＳラン女子高生
ち○ぴょろすぽ○ん
by ドＳなＳラン女子高生
なーんでよそのスレに迷惑をかけに行くのかなん
やっぱりＦラン脳のゆえんなのかしらねんのねん
by ドＳなＳラン女子高生
これやってみ、ド底辺頭脳のオマエには無理だろうな。
24時間待ってやるぞ。
ドツボ13は100発0中
ドツボ14は10発0中
ドツボ15は1発0中
とする。
各々10000発撃ったときドツボの命中数の期待値はいくらか？
感度・特異度を確率変数としてベイズ階層モデルを組めば計算できること投稿しておいた。
最頻値と標準偏差からベータ分布のパラメータを計算するのに連立方程式を解く方が面倒だった。
（キモワロス）
ド底辺仲間のコンプ薬屋はまだですか？
by ドＳなＳラン女子高生
>>88
これやってみ、ド底辺頭脳のオマエには無理だろうな。
24時間待ってやるぞ。
ドツボ13は100発0中
ドツボ14は10発0中
ドツボ15は1発0中
とする。
各々10000発撃ったときドツボの命中数の期待値はいくらか？
国試に合格できなかった還暦ジジイって本当ですかぁ？
アチコチで反日レスしてるみたいですねー
by ドＳなＳラン女子高生
ドツボの宿題まだぁ？
Ｆラン高の宿題はＦランのアナタ自身でされてくださいね
Ｓランのアタシはアナタと違ってヒマではございませんのでＷ
by ドＳなＳラン女子高生
ドツボ問題もできないドアホ。
これやってみ、ド底辺頭脳のオマエには無理だろうな。
24時間待ってやるぞ。
ドツボ13は100発0中
ドツボ14は10発0中
ドツボ15は1発0中
とする。
各々10000発撃ったときドツボの命中数の期待値はいくらか？
この集団の有病率はstanでMCMCして計算すると
mean lower upper
0.7558614 0.5214634 0.9999897
感度60% 特異度90%に固定すると
> pn2ip(95,47)
infected prevalence
75.0000000 0.7894737
病院で管理する医療用麻薬の数量について道に虚偽の届け出をしたとして、
札幌区検は１９日までに麻薬取締法違反の罪で、札幌ひばりが丘病院（札幌
市厚別区）に勤務していた３０代の薬剤師の男と、同病院を運営する医療法人
潤和会（同区）を略式起訴した。札幌簡裁がそれぞれ２０万円以下の罰金刑
を言い渡す見通し。略式命令請求によると、男は２０１５年１１月、道に対し、
院内で使用した麻薬の使用量を偽って届け出るなどしたとされる。潤和会には
法人も罰する同法の両罰規定を適用した。
出典：北海道新聞 平成30年10月19日付
日本人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)
}
全体で何人の感染者がいるかの期待値と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 - (ｍ人全員非感染の確率)
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人')
n(=10)人の中から無作為にm(=2)人選んだらその中に少なくとも一人の感染者がいた。
全体で何人の感染者がいるかの期待値を求めよ。
"
f <- function(p){
# p:感染確率
p1=2*p*(1-p) # 一人だけ感染確率
p2=p^2 # 二人とも感染確率
(1*p1+2*p2)/(p1+p2) # 感染人数の期待値
}
curve(f(x),bty='l',lwd=2,xlab='p')
f(2/3)
# 有名人が感染
" 長時間必要なのて実行非推奨,メモリ16G以下はクラッシュする
library(gmp)
m=18200 # 有名人の数(桜を見る会参加人数)
n=1.268e5 # 日本の人口
x=0:n # 感染者数:x, 非感染数:n-x
pmf=try(1- chooseZ(n-x,m)/chooseZ(n,m)) # 1 - (ｍ人全員非感染の確率)
pdf=try(pmf/sum(pmf)) # 確率密度関数化して
(E=try(sum(x*pdf))) # 期待値を計算
as.numeric(E) # E=63463.27 (m=1000) , E=1154070201/18202=63403.48(m=1.268e5)
"
fn(100,20,print=T)
data.frame(有名人=1:10,期待値=sapply(1:10,function(x) fn(100,x)$mean))
data.frame(有名人=1:10*10,期待値=sapply(1:10*10,function(x) fn(100,x)$mean))
# Pr[B|Ax] n人中x人感染者がいるときにm人中k人感染している確率
library(gmp)
x=0:n
if(quick){ Ev=sum(x*chooseZ(x,k))/sum(chooseZ(x,k))
}else{
pmf=chooseZ(x,k)/chooseZ(n,m)
pdf=pmf/sum(pmf)
Ev=sum(x*pdf)}
# Σ(x*(xCk/nCm))/Σ(xCk/nCm) = Σ(x*(xCk))/Σ(xCk)
if(print) print(Ev)
E=as.numeric(Ev)
return(E)
}
fn4=Vectorize(fn4)
fn4(100,10,0:10,quick=F)
今思うのは反日野党が政権をとってたら
世界はモノスゴクひどいことになっているであろうこと
by 都内Ｓラン女子高生
事故当日に原子力緊急事態宣言を発令、命をかけて原発に突入して手動ベントを約束させる、
自衛隊１０万に派遣決定、
全原発停止、
東電の撤退阻止
と、戦力の逐次投入でなくてもてる力をすべて出し切って日本を守ろうとした。
原発事故の最悪シナリオも専門家に依頼して作成して、それをパワーポイントでレクチャーさせて閣僚で危機感を共有していたという。
一方、安倍は菅直人が海水注入を止めたとデマを流しつづけた。名誉毀損裁判途中でデマメルマガをこっそりと削除。
他スレでお医者さんごっこだぁW
今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの
by 都内Ｓラン女子高生
> summary(js$n.pos.sum)
Min. 1st Qu. Median Mean 3rd Qu. Max.
3888 5094 5951 31450 8608 421156
> hdi(js$n.pos.sum)[1:2]
lower upper
3888 172503
今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの
by 都内Ｓラン女子高生
あぁあ Ｆランジジイだぁ
コテハンを命名してあげましょう
コンプにならってコロナ薬屋とでも名乗るが良いですよ
by 都内Ｓラン女子高生
他スレでお医者さんごっこだぁW
今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの
by 都内Ｓラン女子高生
Fラン統計ジジイ あらため コロナ薬屋が
他スレでお医者さんごっこだぁW
今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの
by 都内Ｓラン女子高生
のデータを使って
P=14895/153581 # 2020/05/04
PCR検査の感度30-70％のモデルM1と感度70-90％のモデルM2のどちらが信憑性があるか、ベイズファクターで計算してみる。
M1は感度が最頻値60%標準偏差10%、M2は最頻値80%標準偏差10%のベータ分布に設定
特異度はいずれも最頻値95%標準偏差2.5%に設定し、有病率は一様分布を仮定
陽性数は、陽性率(P)＝真陽性率＋偽陽性率=有病率=有病率*感度 + (1-有病率)*(1-特異度)の二項分布に従うとする。
事後確率分布は
://i.imgur.com/81bK4KE.png
陽性率P=14895/153581=0.09698465での事後確率分布の密度比(Savage-Dickey density ratio)でベイズファクターを出すと
> d1/d2 # Savage-Dickey densiti ratio = BF12
[1] 1.007722
まあ、ちょっぴり、感度30-70％のモデルの方がいいかも、という結果。
陽性数/検査数の時系列データでもあればもう少し差がでるかもしれん。
PCR検査の感度は30-70%特異度は95-100%として、この患者が非感染者である確率の期待値と信頼区間を計算せよ。
(P.AA_COVID=1-72/202)
# https://jamanetwork.com/journals/jama/fullarticle/2765183
sn=Mv2ab(0.5,0.2^2)
sp=Mv2ab(0.95,0.05^2)
data=list(sn=sn,sp=sp, P.AA_COVID=P.AA_COVID)
modelString='
model{
# P(COVID | anosmia and/or ageusia)
P.COVID_AA = P.AA_COVID*P.COVID/
(P.AA_COVID*P.COVID + P.AA_nonCOVID*(1-P.COVID))
# P(PCR- | aosmia and/or ageusia)
P.nPCR = (1-P.COVID_AA)*spc + P.COVID_AA*(1-sen)
# = TN + FN
P.COVID_nPCR= (1-sen)*P.COVID/((1-sen)*P.COVID + spc*(1-P.COVID))
# P(COVID | PCR-) = P(PCR-|COVID)P(COVID)/P(PCR-)
# = P(PCR-|COVID)P(COVID) /
# (P(PCR-|COVID)P(COVID)+P(PCR-|nonCOVID)P(nonCOVID))
# prior
P.AA_nonCOVID ~ dunif(0.056,0.127)
# https://jamanetwork.com/journals/jama/fullarticle/2765183
P.COVID ~ dbeta(1,1) # prevalence
sen ~ dbeta(sn[1],sn[2]) # sensitivity 0.3-0.7
spc ~ dbeta(sp[1],sp[2]) # specificity 0.9-1.0
}
'
Fラン統計ジジイ あらため コロナ薬屋が
他スレでお医者さんごっこだぁW
ちゃんとコテハン付けないとダメでちよぉぉぉ
今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの
by 都内Ｓラン女子高生
ちゃんとコテハン付けるんでしよぉぉぉお
漢字はまだ読めないんでしたっけねえぇ
by 都内Ｓラン女子高生
BF01 <- function(s1,s2,n1,n2){
library(gmp)
bf=chooseZ(n1,s1)*chooseZ(n2,s2)/chooseZ(n1+n2,s1+s2)*(n1+1)*(n2+1)/(n1+n2+1)
BF01=as.numeric(bf)
return(c('in favor of null-hyposisis : H0'=BF01))
}
(bf01=BF01(s1=424,s2=5416,n1=777,n2=9072))
c('in favor of alternative hyposis : H1' = 1/as.numeric(bf01))
Fラン統計ジジイ あらため コロナ薬屋が
他スレでお医者さんごっこだぁW
ちゃんとコテハン付けないとダメでちよぉぉぉ
今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの
by 都内Ｓラン女子高生
および コンプ薬屋さん って
やっぱりバカなんですねぇぇえぇぇぇ！
ちゃんとコテハン付けないとダメでちよぉぉぉ
今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの
by 都内Ｓラン女子高生
漢字は読めないんでしゅかぁぁあぁぁW
ちゃんとコテハン付けないとダメでちよぉぉぉ
今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの
by 都内Ｓラン女子高生
転売生活でもしてらっしゃるのでしょうか？
http://itest.2ch.sc/egg/test/read.cgi/hosp/1588253497/242
レスを遡ると品不足をアッピールして転売を煽るような書き込みをされてましたね
ということは お仲間の Ｆラン統計ジジイあらため コロナ薬屋さんも ご職業は転売ヤーですか？
by 都内Ｓラン女子高生
次はどんな品を転売予定ですか？
それはそうと
今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの
by 都内Ｓラン女子高生
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)
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 都内Ｓラン女子高生
コロナ薬屋さん
コンプ薬屋さんに共感されてますよ
やっぱりお仲間だったんですね
キンモー
by 都内Ｓラン女子高生
コロナ薬屋 コンプ薬屋が
他スレでお医者さんごっこだぁW
ちゃんとコテハン付けないとダメでちよぉぉぉ
今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの
by 都内Ｓラン女子高生
あらため コロナ薬屋さん
あなたと違って大学は旧帝受験予定ですが
それはそうと昨夜は何を転売するため並んでたんですか？
野糞しながらぼったくり転売は恥ずかしくないのですか？
by 都内Ｓラン女子高生
コロナ患者の収容施設にしよう
?は欠測値。
記憶保持率を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 ?
冬子 ? ? ? ? ? ? ? ? ? ?
あらため コロナ薬屋さん
お仕事は転売ヤーで間違いないですね
頭が悪いのに転売のお仕事の為に統計の勉強をされたんですね
by 都内Ｓラン女子高生
正の相関関係がある と
なるほど PCR信者はバカなんですね
by 都内Ｓラン女子高生
Ｆラン統計ジジイ あらため コロナ薬屋 が
他スレでお医者さんごっこだぁW
ちゃんとコテハン付けないとダメでちよぉぉぉ
今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの
by 都内Ｓラン女子高生
私の周りの高３生にはそんなことを支持する者は誰もいないのですが
どんなＦラン高の生徒の話なんでしょうね
受験は生モノ
勉強時間が足りないのや試験日に体調が悪くて不合格になるのは、本人の実力のうちでしょうに
by 都内Ｓラン女子高生
４人目の冬子の確率を面積比例させてみると
https://i.imgur.com/CHLo8FN.png
Ｆラン統計ジジイ あらため コロナ薬屋さんって
反日パヨクなんですかぁぁあ？
サイテーですね
下水のゾウリムシよりもはるかにサイッテーですね
ダンゴムシのごとくゴロゴロ転がりながら
海より深く反省しなさい
by 都内Ｓラン女子高生
わりと、面倒。
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)
{
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')
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()
使っていないと、ggplot2も忘れてしまうなぁ。
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
Ｆラン統計ジジイ あらため コロナ薬屋さん
および コンプ薬屋さんって
裏口容疑者なんですかぁぁあ？
だと思いましたぁぁあぁぁ
早く出頭するんですよ
あなたたちのレスの全てから
頭の悪さを感じますからねえぇぇぇ
by 都内Ｓラン女子高生
デバッグできたので、事前分布をいろいろ変えてみて再現性を確認。
内視鏡休診で纏まった勉強する時間が取れて( ・∀・)ｲｲ!!
##### 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')
先に開業医板を荒らすから
私がここに降臨して差し上げてるんですよ
それすら判らないくらい知能が低いんですかねぇ
女子高生にブチ切れる異常者のジジイって
どんなホラー映画の地獄絵図ですかＷ
by 都内Ｓラン女子高生
half cauchy にしても一様分布にしても事後分布はほとんど変化はないな。
モデルが優れているからだろう。
お家に帰ってお勉強しましょう。
あ〜あ、Ｆラン統計ジジイは
やっぱりくるくるぱーだったんですねぇ
by 都内Ｓラン女子高生
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))
https://i.imgur.com/TnXbBKY.png
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)
by 都内Ｓラン女子高生
https://asahi.2ch.sc/test/read.cgi/newsplus/1589284722/
東大の児玉名誉教授も同じ意見だな。
https://www.covid19-yamanaka.com/cont4/14.html
色を変えて散布図を作成
https://i.imgur.com/l1eyD6o.png
俺は６割の休業補償だが、銭ゲバ開業医はパート医を首にしているんだろうなぁ。
Ｆラン統計ジジイ あらため コロナ薬屋 が
発狂して連投だぁぁぁW
ちゃんとコテハン付けないとダメでちよぉぉぉ
今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの
by 都内Ｓラン女子高生
内視鏡バイトは休診なので、６割の休業補償はあるけど。
あんたのところは、感染拡大阻止に休診にして職員には休業補償してんの？
まとまった時間がとれるので、https://bayesmodels.com/の英文テキストで独学中。
行き詰まったら数学版で相談。
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
グラフにしてみた。
https://i.imgur.com/OyTi5Jb.png
必要なら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変換の逆関数で出しているな。
sample関数を使ってのbootstrap
https://i.imgur.com/gAIJTEQ.png
"名","本人","家族含む"
"鳩山由紀夫首相",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)
パッケージ 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
アベノミクスでの経済成長率は震災があった民主の頃より低いという現実を認められないんだなぁ。
日本を破壊したいなら、安倍一択。
安倍の長期政権で中国も韓国も笑いが止まらないはず、
日本が衰退する政策ばかりを選択しているから。
野党が反感をもたれて安倍信者が増えることこそ日本破壊が進んで
パヨクの思う壺
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
#--- 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)
"
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)
Ｆラン統計ジジイ あらため コロナ薬屋 が
発狂して連投だぁぁぁW
ちゃんとコテハン付けないとダメでちよぉぉぉ
今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの
by 都内Ｓラン女子高生
data=dataList,n.chains=4,sample=10000,burnin=4000)
coda::gelman.plot(runJags)
codaSamples = as.mcmc.list(runJags)
Ｆラン統計ジジイ あらため コロナ薬屋 が
発狂して連投だぁぁぁW
相変わらず臨床経験や社会経験がゼロなのが
丸わかりのレスですねぇぇぇぇぇ
ちゃんとコテハン付けないとダメでちよぉぉぉ
今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの
by 都内Ｓラン女子高生
id:3+Rla4zY
くるくるぱーが反日野党に肩入れして
自滅するのは自業自得だけど
善良な日本人は巻き込まれたくないですねえ
by 都内Ｓラン女子高生
Ｆラン統計ジジイ あらため コロナ薬屋 が
発狂して連投だぁぁぁW
相変わらず臨床経験や社会経験がゼロなのが
丸わかりのレスですねぇぇぇぇぇ
今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの
by 都内Ｓラン女子高生
現在いる実験室で、様々細胞と遺伝子変異させた細胞に抗がん剤をかけて死亡率をみてます。
抗がん剤の量も比較してるのです
1: 細胞aとbに対して同量の抗がん剤を使用し死亡率を見る場合
2: 同種の細胞に違う濃度の抗がん剤をかけて死亡率を比較する場合。
両者とも何検体かとって平均の比較をする場合、対応するt検定ですか？、それともf値をみた上で対応しないt検定になるのですか？
比率の検定ならΧ（カイ）二乗検定。
Ｆラン統計ジジイ あらため コロナ薬屋 が
発狂して連投だぁぁぁW
身の程を弁えない謎の上から目線で政治談批評
老害アルアルですねぇぇぇぇぇ
今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの
by 都内Ｓラン女子高生
観光立国という亡国政策で衰退途上国。
観光でしか生きていけない国になったときに中国からの出国禁止すれば日本は中国に平伏すしかなくなる。
対数変換して検定した方がいい場合もある。
やっぱりマスコミはクズですね
問題のあった新聞社の社長は直ちに謝罪して
新聞社を解体するべきですね
反日野党は日本の足を引っ張ることしかしてませんね
存在意義がないとしか言えません
だって反日なんだもの
by 都内Ｓラン女子高生
対数変換した場合、検定方法は何使えばいいですか？
平均の検定ならt検定
接客したキャバ嬢Cが人物D発症の２日後に発症していたことがわかった。
キャバ嬢Cは人物Dから移されたと主張して１億円の賠償を求めている。
潜伏期間には幅がありキャバ嬢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')
t検定は二つのことを比べることしかできないですよね？
多重検定やりたいならTukeyとか色々ある。
図の点線が対照のmockで実線が遺伝子変化させた細胞、
横軸が加えた薬剤の濃度です。
ルシRーゼの比を見ています
同じ濃度で、モックと遺伝子変化させた群を比較させたい時は統計は何使えばいいですか？
https://i.imgur.com/VXfJ8TD.jpg
厳しいbonferroni補正したpairwise.t.testとか
Rだとprop.trend.testでできたはず！
すいませんjmpかエクセル統計で計算できるのだと嬉しいです。
kruskal.testで多重比較で有意差確認してから
平均比較ならpairwise.t.test
比率比較ならpairwise.prop.test
補正法は指定できる。
エクセルのマクロは売られている。
使ったこともないけど。
https://bellcurve.jp/ex/function/kruskal.html
Ｆラン統計ジジイ あらため コロナ薬屋 が
発狂して連投だぁぁぁW
それにしても反日野党や反日マスコミは日本の足を引っ張ることしかしてませんね
存在意義がないとしか言えません
だって反日なんだもの
by 都内Ｓラン女子高生
>勤務医が院長のおれの悪口を妻に話してたそうだ
>腹が立ったからクビにする
を登場人物としてみた。
(問題)
新型コロナに勤務医が罹患。
勤務医が発症した翌日に院長が発症、その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)
mean(apply(ABCinfected,1,fn1))
fn2 <- function(x){
x['Cinfected'] < x['Ainfected']
& x['Ainfected'] < x['Binfected']
}
mean(apply(ABCinfected,1,fn2))
( ： )
（ ゜∀゜）ﾉ彡
<( )
ﾉωヽ
自分の尻を両手でバンバン叩きながら白目をむき
从
Дﾟ ） て
( ヾ) )ヾ て
< <
人__人__人__人__人__人__人__人__人__人__人
Σ て
Σ びっくりするほどユートピア！ て人__人_
Σ びっくりするほどユートピア！ て
⌒Ｙ⌒Ｙ⌒Ｙ) て
Ｙ⌒Ｙ⌒Ｙ⌒Ｙ⌒Ｙ⌒Ｙ⌒Ｙ⌒Ｙ⌒Ｙ⌒Ｙ⌒Ｙ⌒
＿_＿＿＿＿＿
|＿＿ ヽ(゜∀゜）ﾉ
|\＿〃´￣￣ ヽ..ヘ( )ミ
| |\,.-〜´￣￣ ω > （∀゜ ）ﾉ
\|∫\ _,. - ､_,. - ､ \ ( ﾍ)
\ \＿＿＿＿＿＿ _\<
\ ||￣￣￣￣￣￣￣ |
\||＿＿＿＿＿＿＿ |
これを10分程続けると妙な脱力感に襲われ、解脱気分に浸れる。
その間隔はt1,t2,..,tn-1としたときにxiが感染源であった確率を算出するプログラムを作れ。
早くアメリカと足並み揃えて中国に損害賠償請求して国民に赦しを乞うべき。
観光立国という亡国政策で日本は衰退途上国。
次の世代の日本人は実習生として中国に出稼ぎがで立場が逆転。
日本人なら過労死自殺しても経営者殺害はしないから酷使される。
観光でしか生きていけない国になったときに中国からの訪日禁止すれば日本は中国に平伏すしかなくなる。
尖閣どころか沖縄を寄越せとか言われても差し出すことになりそうだな。
モデルは簡単なので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)
# 感染日（一人目の発症日：０）
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人発症
Ｆラン統計ジジイ あらため コロナ薬屋は
今日もクルクルパーだよ
忘れないように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];
}
}
'
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))
p
1 0.26
2 0.18
3 0.18
4 0.13
5 0.13
6 0.13
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)
}
# 感染日（一人目の発症日：０）
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))
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]
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
> 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
乱数発生の代わりに数値積分で計算
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)
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)
ゴルゴ問題の答が出せるように勉強したのか？
