臨床統計もおもしろいですよ、その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/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

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

えぇぇえぇぇぇっ

Ｆラン統計ジジイ　あらため　コロナ薬屋さん
および　コンプ薬屋さんって
裏口容疑者なんですかぁぁあ？

だと思いましたぁぁあぁぁ
早く出頭するんですよ
あなたたちのレスの全てから
頭の悪さを感じますからねえぇぇぇ

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/12(火) 08:31:27.05

交絡因子を考慮しての記憶減衰の階層ベイズモデル
デバッグできたので、事前分布をいろいろ変えてみて再現性を確認。
内視鏡休診で纏まった勉強する時間が取れて( ・∀・)ｲｲ!!

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

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

開業医ではないあなたが
先に開業医板を荒らすから
私がここに降臨して差し上げてるんですよ

それすら判らないくらい知能が低いんですかねぇ

女子高生にブチ切れる異常者のジジイって
どんなホラー映画の地獄絵図ですかＷ

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/12(火) 17:39:16.61

>>143
half cauchy にしても一様分布にしても事後分布はほとんど変化はないな。
モデルが優れているからだろう。

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

さて今日も学校は終わり。
お家に帰ってお勉強しましょう。

あ～あ、Ｆラン統計ジジイは
やっぱりくるくるぱーだったんですねぇ

by 都内Ｓラン女子高生

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

# 相関係数の差の検定
n1=150
r1=0.29
n2=120
r2=0.5
pnorm(-abs(atanh(r1)-atanh(r2))/sqrt(1/(n1-3)+1/(n2-3)))

**卵の名無しさん** · 2020/05/13(水) 05:47:04.95

離散量のプロットは重なってしまうと相関がわかりにくくなる。
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))

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

色の濃淡に頻度を反映させるという方法もあるな。

https://i.imgur.com/TnXbBKY.png

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

とりあえず、完成。多分、もっと便利なパッケージがあるんだろうな。

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)

**卵の名無しさん** · 2020/05/13(水) 12:26:32.87

学校の先生に習いましたが、統計的にはPCR検査はむやみにやらないのが正解なんですね

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/13(水) 14:17:48.98

【国民が誤解】安倍晋三「相談・受診の目安＝PCR実施の基準ではない。国民への周知が足りなかった」 [スタス★]
https://asahi.5ch.net/test/read.cgi/newsplus/1589284722/

**卵の名無しさん** · 2020/05/13(水) 14:18:56.03

世界の頭脳は実態把握には大量検査が必要とわかっているなぁ。
東大の児玉名誉教授も同じ意見だな。

https://www.covid19-yamanaka.com/cont4/14.html

**卵の名無しさん** · 2020/05/13(水) 14:55:45.98

>>150
色を変えて散布図を作成

https://i.imgur.com/l1eyD6o.png

**卵の名無しさん** · 2020/05/13(水) 14:57:54.70

検診がメインだったから、内視鏡バイトは休診。
俺は６割の休業補償だが、銭ゲバ開業医はパート医を首にしているんだろうなぁ。

**卵の名無しさん** · 2020/05/13(水) 15:27:10.62

今日も朝から人生が暇しかない
Ｆラン統計ジジイ　あらため　コロナ薬屋　が
発狂して連投だぁぁぁW

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

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

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/13(水) 15:48:20.06

>>156
内視鏡バイトは休診なので、６割の休業補償はあるけど。
あんたのところは、感染拡大阻止に休診にして職員には休業補償してんの？

まとまった時間がとれるので、https://bayesmodels.com/の英文テキストで独学中。
行き詰まったら数学版で相談。

**卵の名無しさん** · 2020/05/13(水) 18:18:49.08

相関係数->FisherのZ変換->ベイズファクター

r=js$r
#The Fisher z-transformation:
r.F <- atanh(r)
fit.posterior <- logspline(r.F)
plot(fit.posterior)
posterior <- dlogspline(atanh(0), fit.posterior)
r.F.prior <- atanh(runif(n.iter,-1,1))
fit.prior <- logspline(r.F.prior)
prior <- dlogspline(atanh(0), fit.prior)
BF10 <- prior/posterior
BF10

**卵の名無しさん** · 2020/05/14(木) 09:44:43.41

まあ、ベイズファクターの計算結果はさほど変わらんな。

> 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

**卵の名無しさん** · 2020/05/14(木) 09:45:47.27

>>159
グラフにしてみた。

https://i.imgur.com/OyTi5Jb.png

**卵の名無しさん** · 2020/05/14(木) 10:03:25.31

あんまり相関係数の信頼区間など気にしていなかった。
必要なら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変換の逆関数で出しているな。

**卵の名無しさん** · 2020/05/14(木) 10:39:23.42

>>161
sample関数を使ってのbootstrap

https://i.imgur.com/gAIJTEQ.png

**卵の名無しさん** · 2020/05/14(木) 11:21:55.13

cat('
"名","本人","家族含む"
"鳩山由紀夫首相",144269,144269
"菅直人副総理兼国家戦略担当相",905,2232
"原口一博総務相",914,1220
"千葉景子法相",3523,3523
"岡田克也外相",3273,8641
"藤井裕久財務相",14356,20214
"川端達夫文部科学相",4024,5583
"長妻昭厚生労働相",0,891
"赤松広隆農相",4864,5934
"直嶋正行経済産業相",3333,3333
"前原誠司国土交通相",741,1441
"小沢鋭仁環境相",2089,4014
"北沢俊美防衛相",309,609
"平野博文官房長官",1195,1875
"中井洽国家公安委員長",1296,1296
"亀井静香金融・郵政担当相",9427,18745
"福島瑞穂消費者・少子化担当相",12734,25000
"仙谷由人行政刷新担当相",1968,3987
', file='shisan.csv')

X=read.csv("shisan.csv")
X
library(boot)
md = boot(X, function(df,idx) median(df[idx,2]), R=1e4)
plot(md)
mn = boot(X, function(df,idx) mean(df[idx,2]), R=1e4)
plot(md)
boxplot(list(median=md$t,mean=mn$t),horizontal = TRUE, las=1)

**卵の名無しさん** · 2020/05/14(木) 11:33:51.60

>>162
パッケージ bootをつかうと

library(boot)
df=data.frame(x,k)
bt=boot(df,function(df,idx) cor(df[idx,1],df[idx,2]),R=1e5)
plot(bt)
HDInterval::hdi(bt$t)[1:2]
BEST::plotPost(as.vector(bt$t),xlab=bquote(cor),compVal = 0)

> HDInterval::hdi(bt$t)[1:2]
[1] -0.06736828 0.31237556

https://i.imgur.com/vCeo1Yn.png

**卵の名無しさん** · 2020/05/14(木) 11:37:20.35

>>156
アベノミクスでの経済成長率は震災があった民主の頃より低いという現実を認められないんだなぁ。

日本を破壊したいなら、安倍一択。
安倍の長期政権で中国も韓国も笑いが止まらないはず、
日本が衰退する政策ばかりを選択しているから。

野党が反感をもたれて安倍信者が増えることこそ日本破壊が進んで
パヨクの思う壺

**卵の名無しさん** · 2020/05/17(日) 15:29:50.58

結局、これが核心部分

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

**卵の名無しさん** · 2020/05/18(月) 07:52:43.74

# https://www.nejm.org/doi/full/10.1056/NEJMoa2001316
#--- incubation period ---
# from Li et al NEJM 2020
# lognormal mean = 5.2
ln.par1 = 1.434065
ln.par2 = 0.6612

x=rlnorm(1e6,ln.par1,ln.par2)
BEST::plotPost(x)
curve(dlnorm(x,ln.par1,ln.par2),add=T)
quantile(x,c(0.025,0.5,0.975))

"
c.onset = c.infected + c.incubation
d.onset = d.infected + d.incubation
c.onset - d.onset = onset.delay
c.infected - d.infected + c.incubation - d.incubation = delay
c.infected - d.infected = onset.delay + d.incubation - c.incubation
"
onset.delay = 2
d.incubation = rlnorm(1e6,ln.par1,ln.par2)
c.incubation = rlnorm(1e6,ln.par1,ln.par2)
infection.delay = onset.delay + d.incubation - c.incubation
BEST::plotPost(infection.delay,compVal = 0)

**卵の名無しさん** · 2020/05/18(月) 13:03:08.05

stancode=
"
data{
real onset_delay;
real ln_par1;
real ln_par2;
}
parameters{
real <lower=0> d_incubation;
real <lower=0> c_incubation;
}
transformed parameters{
real infection_delay = onset_delay + d_incubation - c_incubation;
}
model{
d_incubation ~ lognormal(ln_par1,ln_par2);
c_incubation ~ lognormal(ln_par1,ln_par2);
}

"
model=stan_model(model_code = stancode)
fn.stan <- function(delay){
dataList=list(onset_delay=delay,ln_par1=ln.par1,ln_par2=ln.par2)
fit=sampling(model,data=dataList)
ms=rstan::extract(fit)
mean(ms$infection_delay < 0)
}
fn.stan(2)
onset_delays=0:20
y=sapply(onset_delays,fn.stan)
plot(onset_delays,y, ylab='Pr[ Infected Later ])',axes=F) ; axis(1)

**卵の名無しさん** · 2020/05/18(月) 22:00:40.40

今日も朝から人生が暇しかない
Ｆラン統計ジジイ　あらため　コロナ薬屋　が
発狂して連投だぁぁぁW

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

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

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/19(火) 07:38:54.76

runJags=run.jags('TEMPmodel.txt',monitor=c('p1','p2','diff'),
data=dataList,n.chains=4,sample=10000,burnin=4000)
coda::gelman.plot(runJags)
codaSamples = as.mcmc.list(runJags)

**卵の名無しさん** · 2020/05/19(火) 07:59:02.34

今日も朝から人生が暇しかない
Ｆラン統計ジジイ　あらため　コロナ薬屋　が
発狂して連投だぁぁぁW

相変わらず臨床経験や社会経験がゼロなのが
丸わかりのレスですねぇぇぇぇぇ

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

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

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/21(木) 12:57:01.59

今日のくるくるぱーのIDだああぁぁあ
id:3+Rla4zY

くるくるぱーが反日野党に肩入れして
自滅するのは自業自得だけど
善良な日本人は巻き込まれたくないですねえ
by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/21(木) 17:33:16.75

今日も朝から人生が暇しかない
Ｆラン統計ジジイ　あらため　コロナ薬屋　が
発狂して連投だぁぁぁW

相変わらず臨床経験や社会経験がゼロなのが
丸わかりのレスですねぇぇぇぇぇ

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

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/22(金) 17:43:59.44

すいません、すごい基本的なことなのですが教えてください。

現在いる実験室で、様々細胞と遺伝子変異させた細胞に抗がん剤をかけて死亡率をみてます。
抗がん剤の量も比較してるのです
1: 細胞aとbに対して同量の抗がん剤を使用し死亡率を見る場合
2: 同種の細胞に違う濃度の抗がん剤をかけて死亡率を比較する場合。
両者とも何検体かとって平均の比較をする場合、対応するt検定ですか？、それともf値をみた上で対応しないt検定になるのですか？

**卵の名無しさん** · 2020/05/22(金) 19:39:38.80

>>174
比率の検定ならΧ（カイ）二乗検定。

**卵の名無しさん** · 2020/05/23(土) 07:42:32.91

今日も朝から人生が暇しかない
Ｆラン統計ジジイ　あらため　コロナ薬屋　が
発狂して連投だぁぁぁW

身の程を弁えない謎の上から目線で政治談批評
老害アルアルですねぇぇぇぇぇ

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

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/23(土) 07:57:58.92

日本衰退を願う勢力は安倍一択。
観光立国という亡国政策で衰退途上国。
観光でしか生きていけない国になったときに中国からの出国禁止すれば日本は中国に平伏すしかなくなる。

**卵の名無しさん** · 2020/05/23(土) 07:58:44.46

>>174
対数変換して検定した方がいい場合もある。

**卵の名無しさん** · 2020/05/23(土) 22:14:18.97

テレビを見てるんですが
やっぱりマスコミはクズですね
問題のあった新聞社の社長は直ちに謝罪して
新聞社を解体するべきですね

反日野党は日本の足を引っ張ることしかしてませんね
存在意義がないとしか言えません
だって反日なんだもの

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/24(日) 15:16:08.53

>>178
対数変換した場合、検定方法は何使えばいいですか？

**卵の名無しさん** · 2020/05/24(日) 18:15:55.34

>>180
平均の検定ならt検定

**卵の名無しさん** · 2020/05/25(月) 12:16:59.94

ある人物Dが新型コロナ肺炎に罹患したとする。行動調査によって発症前にキャバクラに行っており
接客したキャバ嬢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')

**卵の名無しさん** · 2020/05/25(月) 12:21:54.11

>>181
t検定は二つのことを比べることしかできないですよね？

**卵の名無しさん** · 2020/05/25(月) 16:55:34.95

>>183
多重検定やりたいならTukeyとか色々ある。

**卵の名無しさん** · 2020/05/25(月) 17:16:47.59

>>184
図の点線が対照のmockで実線が遺伝子変化させた細胞、
横軸が加えた薬剤の濃度です。
ルシフェラーゼの比を見ています
同じ濃度で、モックと遺伝子変化させた群を比較させたい時は統計は何使えばいいですか？
https://i.imgur.com/VXfJ8TD.jpg

**卵の名無しさん** · 2020/05/25(月) 17:49:48.35

好きなの使えば。
厳しいbonferroni補正したpairwise.t.testとか

**卵の名無しさん** · 2020/05/25(月) 18:00:37.69

回帰係数の比較ならコクランアーミテッジだったかな。
Rだとprop.trend.testでできたはず！

**卵の名無しさん** · 2020/05/25(月) 19:07:25.19

皆さんありがとうございます
すいませんjmpかエクセル統計で計算できるのだと嬉しいです。

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

やったね　たえちゃん

**卵の名無しさん** · 2020/05/25(月) 21:59:58.02

Rだと
kruskal.testで多重比較で有意差確認してから
平均比較ならpairwise.t.test
比率比較ならpairwise.prop.test
補正法は指定できる。

エクセルのマクロは売られている。
使ったこともないけど。
https://bellcurve.jp/ex/function/kruskal.html

**卵の名無しさん** · 2020/05/26(火) 07:48:20.04

今日も朝から人生が暇しかない
Ｆラン統計ジジイ　あらため　コロナ薬屋　が
発狂して連投だぁぁぁW

それにしても反日野党や反日マスコミは日本の足を引っ張ることしかしてませんね
存在意義がないとしか言えません
だって反日なんだもの

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/26(火) 12:50:37.25

"
>勤務医が院長のおれの悪口を妻に話してたそうだ
>腹が立ったからクビにする
を登場人物としてみた。
(問題)
新型コロナに勤務医が罹患。
勤務医が発症した翌日に院長が発症、その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)

**卵の名無しさん** · 2020/05/26(火) 12:50:43.88

fn1 <- function(x) min(x)==x['Cinfected']
mean(apply(ABCinfected,1,fn1))

fn2 <- function(x){
x['Cinfected'] < x['Ainfected']
& x['Ainfected'] < x['Binfected']
}
mean(apply(ABCinfected,1,fn2))

**卵の名無しさん** · 2020/05/26(火) 23:37:54.50

まず全裸になり
　　　　　　　　　　　　　( 　：　)
　　　　　　　（゜∀゜）ﾉ彡
　　　　　　 <(　　　)
　　　　　　　　ﾉωヽ

　自分の尻を両手でバンバン叩きながら白目をむき
　　　　　　　　　　　从
　　　　　　　Дﾟ　）　　て
　　　　　　　 (　ヾ) )ヾ　て
　　　　　　　 < <

　　　　　人__人__人__人__人__人__人__人__人__人__人
　　　　Σ　　　　　　　　　　　　　　　　　　　　　　　　　て
　　　　Σ　　びっくりするほどユートピア！　　　　　　　て人__人_
　　　　Σ　　　　　　　　　びっくりするほどユートピア！　　　　　　て
　　　　　⌒Ｙ⌒Ｙ⌒Ｙ)　　　　　　　　　　　　　　　　　　　　　　　　　て
　　　　　　　　　　　　Ｙ⌒Ｙ⌒Ｙ⌒Ｙ⌒Ｙ⌒Ｙ⌒Ｙ⌒Ｙ⌒Ｙ⌒Ｙ⌒Ｙ⌒
　＿_＿＿＿＿＿
　|＿＿　　　　　　ヽ(゜∀゜）ﾉ
　|\＿〃´￣￣ヽ..ヘ(　　 )ミ
　|　|\,.-～´￣￣　　 ω　>　（∀゜）ﾉ
　\|∫\　　 _,. - ､_,. - ､ \　(　　ﾍ)
　　 \ 　 \＿＿＿＿＿＿ _\<
　　　　\ 　||￣￣￣￣￣￣￣ |
　　　　　\||＿＿＿＿＿＿＿ |

これを10分程続けると妙な脱力感に襲われ、解脱気分に浸れる。

**卵の名無しさん** · 2020/05/27(水) 04:43:18.93

x1,x2,...,xnの順に発症
その間隔はt1,t2,..,tn-1としたときにxiが感染源であった確率を算出するプログラムを作れ。

**卵の名無しさん** · 2020/05/27(水) 05:15:26.64

安倍は新コロナは中国発と認めたから春節ウェルカムは外患誘致。
早くアメリカと足並み揃えて中国に損害賠償請求して国民に赦しを乞うべき。

観光立国という亡国政策で日本は衰退途上国。
次の世代の日本人は実習生として中国に出稼ぎがで立場が逆転。
日本人なら過労死自殺しても経営者殺害はしないから酷使される。

観光でしか生きていけない国になったときに中国からの訪日禁止すれば日本は中国に平伏すしかなくなる。
尖閣どころか沖縄を寄越せとか言われても差し出すことになりそうだな。

**卵の名無しさん** · 2020/05/27(水) 06:40:49.96

>>195
モデルは簡単なのでStanやJagsを使わずに完成

WhoInfectedFirst <- function( # 発症間隔から感染源の確率を計算
t=c(1,2), # 発症間隔
k=1e5 # 乱数発生数
){
ln_par1 = 1.434065 # 潜伏期間対数正規分布パラメータ
ln_par2 = 0.6612
n=length(t)+1 # 発症人数

# 潜伏期間
x_incub = matrix(rep(NA,n*k),ncol=n)
for(i in 1:n) x_incub[,i] = rlnorm(k,ln_par1,ln_par2)

# 感染日（一人目の発症日：０）
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人発症