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

**卵の名無しさん** · 2020/03/05(木) 20:17:05.03

　
　内科認定医受験の最低限の知識、
　製薬会社の示してくる臨床データ、
　論文の考察、
　論文を書くときの正当性、
　というのが、臨床統計の今までの目的の大きい部分でしたが、
　
　AI=機械学習の基本も、結局は統計学と確率に支配されます。
　そういう雑多な話をするスレです。
　
※前スレ
臨床統計もおもしろいですよ、その１
https://egg.2ch.net/test/read.cgi/hosp/1493809494/
臨床統計もおもしろいですよ、その2
https://egg.5ch.net/test/read.cgi/hosp/1540905566/

**卵の名無しさん** · 2020/03/15(日) 18:08:53.16

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

**卵の名無しさん** · 2020/03/15(日) 18:09:31.68

シミュレーション解の方が回数上限なしでプログラムできるな。
# 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)

**卵の名無しさん** · 2020/03/15(日) 18:33:54.47

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)

**卵の名無しさん** · 2020/03/15(日) 18:40:27.27

>>41
ド底辺頭脳でも数くらい数えられるんだろ？

1つのサイコロを10000回投げ、出たの目の積が6の倍数になる確率を分数でもとめよ。

**卵の名無しさん** · 2020/03/15(日) 19:19:43.78

>>45
どうでもいいから早くコンプ薬屋を呼んで来なよ
そんなだからモテないんだお

by 都内高学歴ＪＫ

**卵の名無しさん** · 2020/03/15(日) 19:31:13.08

>>46
数も数えられんの？ガイジ？

**卵の名無しさん** · 2020/03/15(日) 21:05:46.74

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

**卵の名無しさん** · 2020/03/15(日) 21:06:08.12

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)

**卵の名無しさん** · 2020/03/15(日) 21:10:18.08

F-value = 2/(1/sensitivity + 1/PPV )　感度とPPVの調和平均

**卵の名無しさん** · 2020/03/16(月) 00:14:20.50

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

**卵の名無しさん** · 2020/03/16(月) 06:35:35.55

>>51
早くＦラン仲間のコンプ薬屋を呼んで来なよ
お勉強してくれるって言ってるお

by 都内高学歴ＪＫ

**卵の名無しさん** · 2020/03/16(月) 06:40:26.18

>>52
ド底辺頭脳でも数くらい数えられるんだろ？

1つのサイコロを10000回投げ、出たの目の積が6の倍数になる確率を分数でもとめよ。

**卵の名無しさん** · 2020/03/16(月) 07:32:37.13

国内で患者数が大幅に増えたときに備えた医療提供体制の確保について
今後、国内で新型コロナウイルス感染症患者数が大幅に増えたときに備え、各都道府県、保健所設置市及び特別区（以下「都道府県等」という。）における外来を受診する患者数等について、
以下の数式を用いて計算いただき、ピーク時の医療需要の目安としてご活用の上、必要な医療提供体制を確保していただくようお願いいたします。

（１）（ピーク時において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

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

**卵の名無しさん** · 2020/03/16(月) 07:32:47.99

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

**卵の名無しさん** · 2020/03/16(月) 07:35:19.73

早くＦラン仲間のコンプ薬屋を呼んで来なよ
童貞仲間でしょぉお

by 都内高学歴ＪＫ

**卵の名無しさん** · 2020/03/16(月) 07:54:48.43

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人の死者予想。

**卵の名無しさん** · 2020/03/16(月) 07:56:56.15

>>56
コンプはこれに答がだせなかったんだが、あんたは出せる？
解答できなきゃ、コンプと同レベルの頭脳と認定してあげよう。
24時間待ってあげる。

問題

某国の新型コロナ感染症の有病率を0.1、PCR検査の感度を0.7　特異度を0.9とする
検査陽性陰性の人を無作為に50人ずつ集めて、カレーを頻回に食べているかを調査した結果が
以下の通りであった。

　　　　　　　　検査陽性　　検査陰性
カレー頻食　a(=18) 　　　　　b(=30)
カレー稀食　c(=32) 　　　　　d(=20)

カレーを頻食すると新型コロナ感染に罹りにくいと結論できるか？　危険率0.05で検定せよ。

**卵の名無しさん** · 2020/03/16(月) 08:12:57.28

国内で患者数が大幅に増えたときに備えた医療提供体制の確保について
今後、国内で新型コロナウイルス感染症患者数が大幅に増えたときに備え、各都道府県、保健所設置市及び特別区（以下「都道府県等」という。）における外来を受診する患者数等について、
以下の数式を用いて計算いただき、ピーク時の医療需要の目安としてご活用の上、必要な医療提供体制を確保していただくようお願いいたします。

（１）（ピーク時において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歳以下の人口は何人と推測されるか？

**卵の名無しさん** · 2020/03/16(月) 13:47:58.64

>>58
評価するには症例数が２桁少ないですわね

それはさておき昔の教授選では積み上げたペーパーの高さを比べたと聞きますが
Ｆランのアナタは採択されたペーパーの１本もあるのかしら？Ｗ

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/03/16(月) 14:14:20.21

やったー　国試無事合格！

**卵の名無しさん** · 2020/03/16(月) 23:11:32.64

>>60
足りないのは標本数でなくてオマエの頭だよ。

**卵の名無しさん** · 2020/03/16(月) 23:31:54.26

数日前、科学技術省の共同予防および制御メカニズムの科学研究グループからの良いニュースがありました。

国家緊急予防および制御薬物工学技術研究センターと深セン第三人民病院LeレイおよびLi英Yチームによって完成された「ファピラビル」。

「新しいコロナウイルス肺炎（COVID-19）患者の安全性と有効性に関する臨床研究」は、予備的な臨床結果に達しています（登録番号：
ChiCTR2000029600）。

研究により、ファビピラビルは、ウイルスクリアランスを促進することにより、新しいコロナウイルス肺炎の進行を緩和することが示されています。研究結果は、中国工学院のジ
ャーナルに提出されました

合計80人の患者がこの臨床試験に登録され、35人の患者がファピラビルで治療され、対照群は年齢、性別、および疾患の重症度が治療群と一致した新冠動脈肺炎の患者45人
でした。 iタブレット治療。薬物投与からウイルス除去までの時間の中央値、治療の14日目の胸部画像の改善率、および安全性を2つのグループ間で比較しました。

結果は、ウイルス除去の観点から、治療後のファピラビル試験群の患者のウイルス除去時間の中央値（ウイルス核酸陰性）は、それぞれ4日と11日であった対照群のそれよりも
有意に短いことを示しました。別の重要な指標として、患者の胸部画像の改善、治療群と対照群の改善率はそれぞれ91.43％と62.22％でした。同時に、対照群と比較し
て、フェラビル治療群は副作用が少なく、忍容性が良好でした。 2つの患者グループのベースライン特性はすべて同等でした。

**卵の名無しさん** · 2020/03/16(月) 23:35:45.65

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

**卵の名無しさん** · 2020/03/16(月) 23:36:27.82

MCMCして信頼区間で比較すると、

https://i.imgur.com/3uZNIwI.png

**卵の名無しさん** · 2020/03/17(火) 00:11:08.71

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

グラフ化なしで統計処理をするな、という例だな。

**卵の名無しさん** · 2020/03/17(火) 00:15:40.58

中央値比較での統計悪用シミュレーションプログラム

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

**卵の名無しさん** · 2020/03/17(火) 00:19:03.10

こんな感じ

https://i.imgur.com/kL3pnn3.png
https://i.imgur.com/3HYkd5O.png

**卵の名無しさん** · 2020/03/17(火) 11:29:17.69

>>64
ド底辺頭脳はサイコロを６回投げたら１の目が１回でているとう風な計算しかできんみたいだな。
信頼区間とか確率分布という概念での思考ができなんよね。
周辺度数固定して二項分布で乱数発生させてχ二乗検定したときのｐ値の分布は以下のようになる。
https://i.imgur.com/5djScD1.png
https://i.imgur.com/RsFr7o4.png
危険率１％で有効といえる確率はほぼ５割であることがわかる。

**卵の名無しさん** · 2020/03/17(火) 11:38:17.66

>>60
エントリーに５以下の数値があるとχ二乗検定だとYatesの補正を使うけど、この程度の標本数は臨床では普通だよ。
中国のアビガンの予備試験も標本数が実薬群ととコントロール群を合わせて80だぞ。
足りないのはオマエの学力だよ。

**卵の名無しさん** · 2020/03/17(火) 12:02:20.58

>>70
まあ　学園カースト最下位のＦランさんが必死ですわね

昨日は先輩方からも国試合格のご報告がたくさんありましたが
Ｆランさんのお知り合いの方たちはいかがですかぁ？Ｗ

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/03/17(火) 12:41:00.33

>>71
答えられない低学力を隠そうと必死だな。
どうやって解くかわからんの？

**卵の名無しさん** · 2020/03/17(火) 19:55:56.01

>>72
わかりました

犯人はアナタですね　Ｆランさん

**卵の名無しさん** · 2020/03/17(火) 21:23:49.50

ようやく完成した

# 球面上の一様分布
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)

**卵の名無しさん** · 2020/03/17(火) 21:46:10.60

# 単位球体内部一様分布の４点で四面体をつくるとき
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)

**卵の名無しさん** · 2020/03/17(火) 21:47:20.28

core i7 Memory 16G はシミュレーションが捗って( ・∀・)ｲｲ!!

**卵の名無しさん** · 2020/03/20(金) 07:04:24.48

臨床医学は確率事象を扱う稼業なので統計処理は必須。

これができない医師はアホといえる。

某国の新型コロナ感染症の有病率を0.1、
PCR検査の感度を0.7　特異度を0.9とする

検査陽性陰性の人を無作為に100人ずつ集めて、カレーを頻回に食べているかを調査した結果が
以下の通りであった。

　　　　　　　　検査陽性　　検査陰性
カレー頻食　　 a(=36) 　　　　b(=60)
カレー稀食　　c(=64) 　　　　d(=40)

カレーを頻食すると新型コロナ感染に罹りにくいと結論できるか？

**卵の名無しさん** · 2020/03/21(土) 10:28:11.17

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

**卵の名無しさん** · 2020/03/21(土) 10:29:01.74

median (interquartile range, IQR)

**卵の名無しさん** · 2020/03/21(土) 17:41:24.46

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

**卵の名無しさん** · 2020/03/22(日) 09:18:45.02

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)

**卵の名無しさん** · 2020/03/22(日) 20:13:20.91

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

**卵の名無しさん** · 2020/03/22(日) 20:14:04.17

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)

**卵の名無しさん** · 2020/03/23(月) 12:03:19.86

何の責任もとれないド底辺ほど
他人の責任にうるさく自分の権利ばかりしつこく要求しがちですわね

コンプ薬屋のことですけど
Fランのアナタも同類ですわよね
早く呼んでくればぁ

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/03/23(月) 13:38:33.35

>>1
よそのスレで騒いでるみたいですけど

Fランのアナタには知性があるとは
医師の方々は考えてはいないようですよ

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/03/23(月) 14:45:03.80

>>1
国試に行き詰まって統計遊びですか

こういうときハゲワロスって言うんでしたっけ？

by ドＳなＳラン女子高生

**卵の名無しさん** · 2020/03/23(月) 14:48:48.41

タバコを吸わないものは自室に禁煙の張り紙など貼らない訳ですが

裏口裏口と盛んにワメキ回るアナタには
身に覚えがあると言うことでよろしいですねＷ

容疑者はオマエだぁっ（コナ○口調で）

by ドＳなＳラン女子高生

**卵の名無しさん** · 2020/03/23(月) 15:00:20.99

アナタの期待値は100発0中だぁっ
ち○ぴょろすぽ○ん

by ドＳなＳラン女子高生

**卵の名無しさん** · 2020/03/23(月) 18:19:32.21

自分で作った統計ごっこのスレがここにあるのに
なーんでよそのスレに迷惑をかけに行くのかなん
やっぱりＦラン脳のゆえんなのかしらねんのねん

by ドＳなＳラン女子高生

**卵の名無しさん** · 2020/03/25(水) 09:16:02.76

>>88
これやってみ、ド底辺頭脳のオマエには無理だろうな。
24時間待ってやるぞ。

ドツボ13は100発0中
ドツボ14は10発0中
ドツボ15は1発0中
とする。
各々10000発撃ったときドツボの命中数の期待値はいくらか？

**卵の名無しさん** · 2020/03/25(水) 09:16:47.51

正確な感度も特異度もわからないから有病率は計算できないという主張が数学板でされたから、
感度・特異度を確率変数としてベイズ階層モデルを組めば計算できること投稿しておいた。
最頻値と標準偏差からベータ分布のパラメータを計算するのに連立方程式を解く方が面倒だった。

**卵の名無しさん** · 2020/03/25(水) 14:59:53.46

きゃー、統計ごっこのクルクルパー来たーっW
（キモワロス）

ド底辺仲間のコンプ薬屋はまだですか？

by ドＳなＳラン女子高生

**卵の名無しさん** · 2020/03/25(水) 15:55:21.10

>>92
>>88
これやってみ、ド底辺頭脳のオマエには無理だろうな。
24時間待ってやるぞ。

ドツボ13は100発0中
ドツボ14は10発0中
ドツボ15は1発0中
とする。
各々10000発撃ったときドツボの命中数の期待値はいくらか？

**卵の名無しさん** · 2020/03/26(木) 10:22:15.80

>>93
国試に合格できなかった還暦ジジイって本当ですかぁ？
アチコチで反日レスしてるみたいですねー

by ドＳなＳラン女子高生

**卵の名無しさん** · 2020/03/26(木) 10:47:49.20

>>94
ドツボの宿題まだぁ？

**卵の名無しさん** · 2020/03/26(木) 11:18:27.82

>>95
Ｆラン高の宿題はＦランのアナタ自身でされてくださいね
Ｓランのアタシはアナタと違ってヒマではございませんのでＷ

by ドＳなＳラン女子高生

**卵の名無しさん** · 2020/03/26(木) 19:27:21.07

>>96
ドツボ問題もできないドアホ。
これやってみ、ド底辺頭脳のオマエには無理だろうな。
24時間待ってやるぞ。

ドツボ13は100発0中
ドツボ14は10発0中
ドツボ15は1発0中
とする。
各々10000発撃ったときドツボの命中数の期待値はいくらか？

**卵の名無しさん** · 2020/03/27(金) 07:08:38.46

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

**卵の名無しさん** · 2020/04/04(土) 05:13:53.27

麻薬取締法違反　札幌の薬剤師と法人を略式起訴

病院で管理する医療用麻薬の数量について道に虚偽の届け出をしたとして、
札幌区検は１９日までに麻薬取締法違反の罪で、札幌ひばりが丘病院（札幌
市厚別区）に勤務していた３０代の薬剤師の男と、同病院を運営する医療法人
潤和会（同区）を略式起訴した。札幌簡裁がそれぞれ２０万円以下の罰金刑
を言い渡す見通し。略式命令請求によると、男は２０１５年１１月、道に対し、
院内で使用した麻薬の使用量を偽って届け出るなどしたとされる。潤和会には
法人も罰する同法の両罰規定を適用した。
出典：北海道新聞　平成30年10月19日付

**卵の名無しさん** · 2020/04/07(火) 08:43:20.38

晋型コロナ肺炎に感度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)
}

**卵の名無しさん** · 2020/04/22(水) 06:24:37.56

"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 - (ｍ人全員非感染の確率)
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人')

**卵の名無しさん** · 2020/04/22(水) 06:25:30.78

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

**卵の名無しさん** · 2020/04/22(水) 06:26:27.83

fn4 <- function(n,m,k,quick=TRUE,print=FALSE){
# 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)

**卵の名無しさん** · 2020/04/24(金) 13:08:57.65

休校続きでお勉強くらいしかすることがございません

今思うのは反日野党が政権をとってたら
世界はモノスゴクひどいことになっているであろうこと

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/04/26(日) 12:36:58.52

今にして思えば菅直人は実に優秀だったな。　　　　　
事故当日に原子力緊急事態宣言を発令、命をかけて原発に突入して手動ベントを約束させる、

自衛隊１０万に派遣決定、
全原発停止、
東電の撤退阻止

と、戦力の逐次投入でなくてもてる力をすべて出し切って日本を守ろうとした。
原発事故の最悪シナリオも専門家に依頼して作成して、それをパワーポイントでレクチャーさせて閣僚で危機感を共有していたという。

一方、安倍は菅直人が海水注入を止めたとデマを流しつづけた。名誉毀損裁判途中でデマメルマガをこっそりと削除。

**卵の名無しさん** · 2020/04/27(月) 08:26:59.50

今日も朝からクルクルパーのFラン統計ジジイが
他スレでお医者さんごっこだぁW

今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/04/30(木) 07:03:38.11

推定陽性者数：
> 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

**卵の名無しさん** · 2020/04/30(木) 23:55:05.47

統計ジジイとコンプ薬屋はＦラン仲間

今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/01(金) 19:26:10.83

いまだに安倍が日本人の敵と気づかないドアホかよ。

**卵の名無しさん** · 2020/05/01(金) 20:58:05.28

>>109
あぁあ　Ｆランジジイだぁ
コテハンを命名してあげましょう

コンプにならってコロナ薬屋とでも名乗るが良いですよ

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/02(土) 11:49:24.12

今日も朝からクルクルパーのFラン統計ジジイが
他スレでお医者さんごっこだぁW

今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/05(火) 12:24:58.63

今日も朝からクルクルパーの
Fラン統計ジジイ　あらため　コロナ薬屋が
他スレでお医者さんごっこだぁW

今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/05(火) 15:39:04.77

https://toyokeizai.net/sp/visual/tko/covid19/
のデータを使って
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％のモデルの方がいいかも、という結果。
陽性数/検査数の時系列データでもあればもう少し差がでるかもしれん。

**卵の名無しさん** · 2020/05/06(水) 13:43:13.77

味覚障害がある患者をPCR検査に回したら陰性であった。
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
}
'

**卵の名無しさん** · 2020/05/06(水) 14:24:40.20

今日も朝からクルクルパーの
Fラン統計ジジイ　あらため　コロナ薬屋が
他スレでお医者さんごっこだぁW

ちゃんとコテハン付けないとダメでちよぉぉぉ

今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/06(水) 14:44:32.79

ボクちゃん
ちゃんとコテハン付けるんでしよぉぉぉお

漢字はまだ読めないんでしたっけねえぇ

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/06(水) 18:04:44.54

# s1,s2,n1,n2 -> BF10
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))

**卵の名無しさん** · 2020/05/08(金) 12:35:41.21

今日も朝から人生が暇しかないの
Fラン統計ジジイ　あらため　コロナ薬屋が
他スレでお医者さんごっこだぁW

ちゃんとコテハン付けないとダメでちよぉぉぉ

今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/08(金) 12:44:01.53

Fラン統計ジジイ　あらため　コロナ薬屋さん
および　コンプ薬屋さん　って
やっぱりバカなんですねぇぇえぇぇぇ！

ちゃんとコテハン付けないとダメでちよぉぉぉ

今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/08(金) 13:28:12.61

Fラン統計ジジイ　あらため　コロナ薬屋さん
漢字は読めないんでしゅかぁぁあぁぁW

ちゃんとコテハン付けないとダメでちよぉぉぉ

今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/08(金) 15:02:48.93

ところで　コンプ薬屋という方は
転売生活でもしてらっしゃるのでしょうか？
http://itest.5ch.net/egg/test/read.cgi/hosp/1588253497/242

レスを遡ると品不足をアッピールして転売を煽るような書き込みをされてましたね

ということは　お仲間の　Ｆラン統計ジジイあらため　コロナ薬屋さんも　ご職業は転売ヤーですか？

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/08(金) 18:33:14.28

それで　転売ヤーの　コロナ薬屋さん　コンプ薬屋さんは

次はどんな品を転売予定ですか？

それはそうと
今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/08(金) 19:48:15.39

library(rjags)
s1=22;s2=10;n1=158;n2=78
c(s1/n1,s2/n2,s1/n1-s2/n2)
shape1=shape2=0.5
data=list(s1=s1,s2=s2,n1=n1,n2=n2,shape1=shape1,shape2=shape2)
modelString='
model{
# data
s1 ~ dbin(p1,n1)
s2 ~ dbin(p2,n2)
# model
p1 ~ dbeta(shape1,shape2)
p2 ~ dbeta(shape1,shape2)
delta = p1 - p2
# priors
pri.delta= pri.p1 - pri.p2
pri.p1 ~ dbeta(shape1,shape2)
pri.p2 ~ dbeta(shape1,shape2)
}
'
writeLines(modelString,'tmp.txt')
n.iter=1e4 ; thin=1
jagsModel=jags.model('tmp.txt',data=data,n.chains=4,n.adapt=n.iter/5)
update(jagsModel)
codaSamples=coda.samples(jagsModel,c('delta','pri.delta'),n.iter,thin)
plot(codaSamples)
js=as.data.frame(as.matrix(codaSamples))
BEST::plotPost(js$delta,compVal = 0)

**卵の名無しさん** · 2020/05/08(金) 19:48:21.41

library(polspline)
fit=logspline(js$delta)
pri.fit=logspline(js$pri.delta)
xlim=c(-1,1) ; ylim=c(0,9)
plot(fit,xlim=xlim, ylim=ylim,bty='l',lty=2)
par(new=T)
plot(pri.fit,xlim=xlim,ylim=ylim)

dpost=dlogspline(0,fit)
dprio=dlogspline(0,pri.fit)
points(0,dpost)
points(0,dprio,pch=19)
dpost/dprio

**卵の名無しさん** · 2020/05/08(金) 20:57:51.65

転売ヤーの　コロナ薬屋さん　コンプ薬屋さん
それが次の転売品ですか？

高音ぼったくりで通販サイトを荒らして
恥ずかしくないのですか？
転売品仕入れの時は夜中から並んで野糞をしてるって本当ですか？

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/08(金) 23:58:03.61

http://itest.5ch.net/egg/test/read.cgi/hosp/1588253497/296
コロナ薬屋さん
コンプ薬屋さんに共感されてますよ

やっぱりお仲間だったんですね
キンモー

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/09(土) 07:25:33.37

今日も朝から人生が暇だけの
コロナ薬屋　コンプ薬屋が
他スレでお医者さんごっこだぁW

ちゃんとコテハン付けないとダメでちよぉぉぉ

今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/09(土) 08:18:45.43

それで　私立薬学部卒の　Ｆラン統計ジジイ
あらため　コロナ薬屋さん

あなたと違って大学は旧帝受験予定ですが

それはそうと昨夜は何を転売するため並んでたんですか？
野糞しながらぼったくり転売は恥ずかしくないのですか？

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/10(日) 00:00:00.36

全てのパチンコ屋を営業停止にして
コロナ患者の収容施設にしよう

**卵の名無しさん** · 2020/05/10(日) 05:08:19.48

3人の女子大生に18種類の体位を教えてtヶ月後に覚えていた体位の数は以下の通りである.
?は欠測値。
記憶保持率をexp(-α*t)+βのモデルに当て嵌めて
4人めの冬子を含めて欠測値を推定せよ。

t=1 2 4 7 12 21 35 59 99 200
治子 18 18 16 13 9 6 4 4 4 ?
奈津子 17 13 9 6 4 4 4 4 4 ?
亜希子 14 10 6 4 4 4 4 4 4 ?
冬子 ? ? ? ? ? ? ? ? ? ?

**卵の名無しさん** · 2020/05/10(日) 08:38:11.59

それで　私立薬学部卒の　Ｆラン統計ジジイ
あらため　コロナ薬屋さん

お仕事は転売ヤーで間違いないですね
頭が悪いのに転売のお仕事の為に統計の勉強をされたんですね

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/10(日) 08:55:40.13

統計上は　PCR件数と　死亡数は
正の相関関係がある　と

なるほど　PCR信者はバカなんですね

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/11(月) 14:17:46.13

今日も朝から人生が暇しかない
Ｆラン統計ジジイ　あらため　コロナ薬屋　が
他スレでお医者さんごっこだぁW

ちゃんとコテハン付けないとダメでちよぉぉぉ

今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/11(月) 15:25:59.72

ネットを見ると、大学の９月入学を支持する高３生の泣き言、みたいなニュースがありますね

私の周りの高３生にはそんなことを支持する者は誰もいないのですが
どんなＦラン高の生徒の話なんでしょうね

受験は生モノ
勉強時間が足りないのや試験日に体調が悪くて不合格になるのは、本人の実力のうちでしょうに

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/11(月) 19:32:29.89

>>130
４人目の冬子の確率を面積比例させてみると
https://i.imgur.com/CHLo8FN.png

**卵の名無しさん** · 2020/05/11(月) 20:28:09.06

えぇぇえぇぇぇっ

Ｆラン統計ジジイ　あらため　コロナ薬屋さんって
反日パヨクなんですかぁぁあ？

サイテーですね
下水のゾウリムシよりもはるかにサイッテーですね
ダンゴムシのごとくゴロゴロ転がりながら
海より深く反省しなさい

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/11(月) 22:43:11.83

>>130
わりと、面倒。
library(rjags)
t <- c(1, 2, 4, 7, 12, 21, 35, 59, 99, 200)
nt <- length(t)
slist <- 1:4
ns <- length(slist)

k <- matrix(c(18, 18, 16, 13, 9, 6, 4, 4, 4, NA,
17, 13, 9, 6, 4, 4, 4, 4, 4, NA,
14, 10, 6, 4, 4, 4, 4, 4, 4, NA,
NA, NA, NA, NA,NA,NA,NA,NA,NA, NA), nrow=ns, ncol=nt, byrow=T)
colnames(k)=paste0('t',as.character(t))
rownames(k)=paste0('subject',1:4)
n <- 18
dataList=list(t=t,k=k,nt=nt,ns=ns,n=n)

**卵の名無しさん** · 2020/05/11(月) 22:43:15.93

cat("model
{
for(i in 1:ns){
for(j in 1:nt){
k[i,j] ~ dbin(theta[i,j],n)
predk[i,j] ~ dbin(theta[i,j],n)
}
}
for(i in 1:ns){
for(j in 1:nt){
theta[i,j] <- min(1 - 1e-16, exp(-alpha[i]*t[j]) + beta[i])
} # min(1, ....) => Node inconsistent with parents -- ERROR
}
for(i in 1:ns){
alpha[i] ~ dbeta(shape1a,shape2a) # dnorm(mu.a,lambda.a)T(0,1)
beta[i] ~ dbeta(shape1b,shape2b) # dnorm(mu.b,lambda.b)T(0,1)
}
shape1a ~ dt(0,pow(2.5,-2),1)T(0,) # mu.a ~ dbeta(1,1)
shape2a ~ dt(0,pow(2.5,-2),1)T(0,) # mu.b ~ dbeta(1,1)
shape1b ~ dt(0,pow(2.5,-2),1)T(0,) # lambda.a ~ dgamma(0.001,0.001)T(0.001,)
shape2b ~ dt(0,pow(2.5,-2),1)T(0,) # lambda.b ~ dgamma(0.001,0.001)T(0.001,)
}
", file='tmp.txt')

**卵の名無しさん** · 2020/05/11(月) 22:44:15.44

pk4=NULL
for(i in 1:10){
pk4=cbind(pk4,eval(str2lang(paste0("js$\'predk[4," ,i, "]\'"))))
}
colnames(pk4)=c(paste0('t0',1:9),'t10')
boxplot(pk4)
head(pk4)
df4=as.data.frame(pk4)
head(df4)
predk4=tidyr::pivot_longer(df4,col=1:10)
head(predk4)
colnames(predk4)=c('time','items')
library(ggplot2)
g <- ggplot(predk4,aes(time,items))
g + geom_violin()
g + geom_boxplot()

**卵の名無しさん** · 2020/05/11(月) 22:45:23.11

reshape2 が　tidyr::pivot_longer(df4,col=1:10)　になっていた。
使っていないと、ggplot2も忘れてしまうなぁ。

**卵の名無しさん** · 2020/05/12(火) 07:14:59.85

g <- ggplot(predk4,aes(time,items))
g + geom_boxplot()
https://i.imgur.com/krliw0d.png

バイオリンプロットというよりナメクジプロットだな。
g + geom_violin()
https://i.imgur.com/fYZHA53.png

histの$densityを使って面積比で反映させてみた
https://i.imgur.com/bIigw6R.png

標準のboxplotでも十分な気がする。
https://i.imgur.com/L210Sf8.png