20/01/28 15:48:52.17 G+B9uPjF.net
## 2x2 tableからsensitivity,specificity,PPV,NPVの信頼区間を計算
cont2ci<- function(TP,FP,FN,TN){ # contingency table 2 confidence interval
library(BEST)
sensitivity=TP/(TP+FN)
specificity=TN/(TN+FP)
pLR=TP/FP # sensitivity/(1-specificity)
nLR=FN/TN # (1-TP)/(1-FP) = (1-sensitivity)/specificity
PPV=TP/(TP+FP)
NPV=TN/(TN+FN)
par(mfrow=c(3,2))
col=sample(colours(),1)
plotPost(sensitivity,showMode = F,col=col)
plotPost(specificity,showMode = F,col=col)
plotPost(PPV,showMode = F,col=col)
plotPost(NPV,showMode = F,col=col)
plotPost(pLR,showMode = F,col=col)
plotPost(nLR,showMode = F,col=col)
sn=hdi(sensitivity)[1:2]
sp=hdi(specificity)[1:2]
ppv=hdi(PPV)[1:2]
npv=hdi(NPV)[1:2]
plr=hdi(pLR)[1:2]
nlr=hdi(nLR)[1:2]
par(mfrow=c(1,1))
rbind(sn,sp,ppv,npv,plr,nlr)
}
900:卵の名無しさん
20/01/28 15:50:33.03 G+B9uPjF.net
# ポアソンモデル
k=1e5
TP=rpois(k,43)
FP=rpois(k,2)
FN=rpois(k,35)
TN=rpois(k,25)
cont2ci(TP,FP,FN,TN)
# 二項分布モデル(周辺度数固定)
k=1e5
TP=rbinom(k,43+35,43/(43+35))
FP=rbinom(k, 2+25, 2/( 2+25))
FN=rbinom(k,35+43,35/(35+43))
TN=rbinom(k,25+ 2,25/(25+ 2))
cont2ci(TP,FP,FN,TN)
# 二項分布モデル(総数固定)
k=1e5
mat=matrix(c(43,2,35,25),ncol=2)
colnames(mat)=c("逆転陽性","逆転陰性")
rownames(mat)=c("除菌成功","除菌失敗")
mat
N=apply(mat,1,sum) # 78 vs 27
N.total=sum(mat)
Succ=rbinom(k,N.total,78/N.total)
Fail=rbinom(k,N.total,27/N.total)
FP=rbinom(k,Fail, 2/Fail)
FN=rbinom(k,Succ,35/Succ)
TN=Fail - FP
TP=Succ - FN
cont2ci(TP,FP,FN,TN)
901:コンブ薬屋
20/01/28 17:55:30 xazq0gTe.net
では次の問題です
この病棟における入院期間の平均値と治癒率を求めよ
URLリンク(togetter.com)
902:卵の名無しさん
20/01/28 19:19:47.40 vtD6D8ed.net
0がエントリーにあると扱いが難しいな。
rm(list=ls())
mat.smoke=cbind(c(0,3),c(13,28)-c(0,3))
mat.smoke # 喫煙のICU入院のクロス表
Epi::twoby2(mat.smoke)
fmsb::oddsratio(mat.smoke)
library(BayesFactor)
library(BEST)
k=1e3
bfp=contingencyTableBF(mat.smoke,'poisson',post=T,iter=k)
matp=as.data.frame(bfp) # ポアソン分布で乱数化
apply(matp,2,mean)
bfi=contingencyTableBF(mat.smoke,'indepMulti','rows',post=T,iter=k)
bfi=as.data.frame(bfi) # 行の周辺度数固定の二項分布で乱数化
mati=bfi[,1:4]*sum(mat.smoke)
apply(mati,2,mean)
bfj=contingencyTableBF(mat.smoke,'jointMulti',post=T,iter=k)
matj=as.matrix(bfj)*sum(mat.smoke) # 総数数固定の二項分布で乱数化
apply(matj,2,mean)
903:卵の名無しさん
20/01/28 19:19:58.78 vtD6D8ed.net
mat2pf <- function(x){
fisher.test(matrix(round(x),nrow=2))$p.value # roundで整数化してFisherテスト
}
mat2pc <- function(x){
chisq.test(matrix(x,nrow=2))$p.value # χ二乗検定
}
ppf=apply(matp,1,mat2pf)
plotPost(ppf,compVal=0.05)
ppc=apply(matp,1,mat2pc)
plotPost(ppc,compVal=0.05)
pif=apply(mati,1,mat2pf)
plotPost(pif,compVal=0.05)
pic=apply(mati,1,mat2pc)
plotPost(pic,compVal=0.05)
pjf=apply(matj,1,mat2pf)
plotPost(pjf,compVal=0.05)
pjc=apply(matj,1,mat2pc)
plotPost(pjc,compVal=0.05)
source('tools.R')
jags4prop(0,3,13,28)
904:卵の名無しさん
20/01/28 19:22:05.70 vtD6D8ed.net
リンクは踏まないから、こういう風に問題を書けよ。
"某女子大には決して嘘をつかない学生と必ず嘘をつく学生がいることがわかっている。
この女子大の学生(嘘つきかどうかは不明)から
「あなたのいうことが正しければ手コキかフェラをしてあげる、間違っていれば何もしてあげない」と言われた。
この女子大生にフェラをしてもらうには何と言えばいいか?
905:卵の名無しさん
20/01/28 19:47:42.30 xazq0gTe.net
ここは下品ないんたねっとですね
まさにこのスレはド底辺と呼ぶにふさわしいですねW
906:卵の名無しさん
20/01/28 20:34:26 vtD6D8ed.net
これで、
# 2x2のクロス表からstanを使ってMCMCサンプリングしてFisherテストしてp.valueの配列からその分布を表示させる
library(rstan)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
contingency2stan2pf <- function(mat,alpha=0.05,print=T){
N=apply(mat,1,sum)
data <- list(mat=mat,N=N)
stanModel=readRDS('contingency.rds')
fit=sampling(stanModel,data=data, chain=1, iter=25000,warmup=5000)
pars=c('y1','n1_y1','y0','n0_y0','d','RR','OR','NNT')
if(print) print(fit, pars=pars,digits=4,prob=c(0.025,0.5,0.975))
ms=rstan::extract(fit)
p.sim=NULL
for(i in 1:length(p11)){
MAT=round(with(ms,matrix(c(y1[i],n1_y1[i],y0[i],n0_y0[i]),ncol=2,byrow=TRUE)))
p.sim=append(p.sim,fisher.test(MAT)$p.value)
}
BEST::plotPost(p.sim,compVal = 0.05, xlab='p.value')
cat('Pr(p.value < ',alpha,') = ', mean(p.sim<alpha),'\n')
cat('Original p.value = ',fisher.test(mat)$p.value,'\n')
print(summary(p.sim))
invisible(p.sim)
}
サンプル数を増やしてpを小さくして信用させるという手法を見破れるかな。
907:卵の名無しさん
20/01/28 20:35:30 vtD6D8ed.net
>>854
ちゃんと答えないとフェラしてもらえないよ。
908:卵の名無しさん
20/01/28 20:38:30 vtD6D8ed.net
> (mat.steroid=cbind(c(4,2),c(9,32)-c(4,2))) # ステロイドと死亡数
[,1] [,2]
[1,] 4 5
[2,] 2 30
> contingency2stan2pf(mat.steroid,p=F)
Pr(p.value < 0.05 ) = 0.7024
Original p.value = 0.01481455782
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0000121509 0.0055484536 0.0148145578 0.0481260381 0.0611450791 1.0000000000
909:卵の名無しさん
20/01/28 21:01:49.02 vtD6D8ed.net
>>855
フィッシャー検定は非負整数しか許さないから、こんな感じで
URLリンク(i.imgur.com)
グラフにすると途切れができてしまうのが美しくない。
χ二乗検定版を作ってみよう。
910:卵の名無しさん
20/01/28 21:03:52.80 IbKEPttD.net
国試浪人の素人童貞のスレはここですか?
911:卵の名無しさん
20/01/28 21:07:29.78 vtD6D8ed.net
>>858
# χ二乗テスト版
contingency2stan2pc <- function(mat,alpha=0.05,print=T){
N=apply(mat,1,sum)
data <- list(mat=mat,N=N)
stanModel=readRDS('contingency.rds')
fit=sampling(stanModel,data=data, chain=1, iter=25000,warmup=5000)
pars=c('y1','n1_y1','y0','n0_y0','d','RR','OR','NNT')
if(print) print(fit, pars=pars,digits=4,prob=c(0.025,0.5,0.975))
ms=rstan::extract(fit)
p.sim=NULL
for(i in 1:length(p11)){
MAT=with(ms,matrix(c(y1[i],n1_y1[i],y0[i],n0_y0[i]),ncol=2,byrow=TRUE))
p.sim=append(p.sim,chisq.test(MAT)$p.value)
}
BEST::plotPost(p.sim,compVal = 0.05, xlab='p.value')
cat('Pr(p.value < ',alpha,') = ', mean(p.sim<alpha),'\n')
cat('Original p.value = ',chisq.test(mat)$p.value,'\n')
print(summary(p.sim))
invisible(p.sim)
}
(mat.smoke=cbind(c(0,3),c(13,28)-c(0,3))) # 喫煙とICU入院のクロス表
contingency2stan2pc(mat.smoke,p=F)
(mat.steroid=cbind(c(4,2),c(9,32)-c(4,2))) # ステロイドと死亡数
p.sim=contingency2stan2pc(mat.steroid,p=F)
sapply(c(0.0001,0.001,0.01,0.05),function(x) mean(p.sim<x))
lines(density(p.sim))
この方が美しいな。
URLリンク(i.imgur.com)
912:卵の名無しさん
20/01/28 21:08:19.29 vtD6D8ed.net
>>859
これできないとフェラしてもらえないよ。
"某女子大には決して嘘をつかない学生と必ず嘘をつく学生がいることがわかっている。
この女子大の学生(嘘つきかどうかは不明)から
「あなたのいうことが正しければ手コキかフェラをしてあげる、間違っていれば何もしてあげない」と言われた。
この女子大生にフェラをしてもらうには何と言えばいいか?
913:卵の名無しさん
20/01/28 21:10:44.10 vtD6D8ed.net
人類の三大発明は
言語・貨幣・フェラチオである。
by 吉田戦車
914:卵の名無しさん
20/01/28 21:52:45.89 cUmejCEO.net
お兄ちゃん、イサキは、イサキは獲れたのっ?
915:卵の名無しさん
20/01/28 21:57:08.65 vtD6D8ed.net
比率の比較でなくて平均値の比較のときはp値の信頼区間算出はブートストラップを使えば簡単にだせるな。
知能指数が平均101と99の各々500人からなる集団の比較
# t検定(等分散)
T.test=function(n1,n2,m1,m2,sd1,sd2){
SE12=sqrt((1/n1+1/n2)*((n1-1)*sd1^2+(n2-1)*sd2^2)/((n1-1)+(n2-1)))
T=(m1-m2)/SE12
p.value=2*pt(abs(T),n1-1+n2-1,lower.tail = FALSE)
return(c(p.value=p.value,t.statistic=T))
}
n=500
a=99
b=101
A=a+scale(rnorm(n))*15
B=b+scale(rnorm(n))*15
T.test(n,n,a,b,15,15)
t.test(A,B,var=T)$p.value
sim <- function(x=n) t.test(sample(A,x,rep=T),sample(B,x,rep=T))$p.value
k=1e4
pv=replicate(k,sim())
summary(pv)
quantile(pv,prob=c(0.005,0.025,0.5,0.975,0.995))
BEST::plotPost(pv,compVal = 0.05,showMode = T,breaks=50)
916:卵の名無しさん
20/01/28 22:00:00.87 vtD6D8ed.net
>>863
釣り人(嘘つきかどうかは不明)から
「あなたのいうことが正しければイサキかイサキの餌のどちらかをあげる、間違っていれば何もあげない」と言われた。
この釣り人からイサキをもらうには何と言えばいいか?
917:卵の名無しさん
20/01/29 05:29:27 qmYv2xxr.net
少し一般化、復元抽出するか否
918:かや個数が不均等でも計算できるように改変。 # 母集団Na,Nb個 平均Ma,Mb 標準偏差SDa,SDbからna,nbを抽出しt検定比較ときのp値の分布 pCI <- function(Na=500,Nb=Na,Ma=99,Mb=101,SDa=15,SDb=SDa,na=50,nb=na, rep=FALSE,k=1e3,print=TRUE){ A=Ma+scale(rnorm(Na))*SDa B=Mb+scale(rnorm(Nb))*SDb sim <- function(m) t.test(sample(A,na,replace=rep),sample(B,nb,replace=rep))$p.value pv=replicate(k,sim()) if(print){ print(T.test(n,n,a,b,15,15)) print(summary(pv)) BEST::plotPost(pv,compVal = 0.05,showMode = T,breaks=50)} invisible(pv) }
919:卵の名無しさん
20/01/29 05:46:52.67 qmYv2xxr.net
デフォルト設定では母集団から1/10の標本を抽出するように変更
# 母集団Na,Nb個 平均Ma,Mb 標準偏差SDa,SDbからna,nbを抽出しt検定比較ときのp値の分布
pCI <- function(Na=500,Nb=Na,Ma=99,Mb=101,SDa=15,SDb=SDa,na=Na/10,nb=Nb/10,
rep=FALSE,k=1e3,print=TRUE){
A=Ma+scale(rnorm(Na))*SDa
B=Mb+scale(rnorm(Nb))*SDb
sim <- function(m) t.test(sample(A,na,replace=rep),sample(B,nb,replace=rep))$p.value
pv=replicate(k,sim())
if(print){
print(t.test(A,B))
print(summary(pv))
BEST::plotPost(pv,compVal = 0.05,showMode = T,breaks=50)}
invisible(pv)
}
920:卵の名無しさん
20/01/29 05:53:26.38 qmYv2xxr.net
標準大学の新入生の知能指数の平均を100とする。
ド底辺シリツ医大の新入生の知能指数の平均が85であった。
各々から1/10を無作為抽出して知能指数をt検定したときのp値の期待値、中央値を求めよ。
また、p値が0.05以上になってド底辺シリツ医の新入生の知能指数は統計的に有意差はないと主張できる確率はいくらか?
# 母集団Na,Nb個 平均Ma,Mb 標準偏差SDa,SDbからna,nbを抽出しt検定比較ときのp値の分布
pCI <- function(Na=500,Nb=Na,Ma=99,Mb=101,SDa=15,SDb=SDa,na=Na/10,nb=Nb/10,
rep=FALSE,k=1e3,print=TRUE){ # IQ:m=100,sd=15
A=Ma+scale(rnorm(Na))*SDa
B=Mb+scale(rnorm(Nb))*SDb
sim <- function(m) t.test(sample(A,na,replace=rep),sample(B,nb,replace=rep))$p.value
pv=replicate(k,sim())
if(print){
print(t.test(A,B))
print(summary(pv))
BEST::plotPost(pv,compVal = 0.05,showMode = T,breaks=50)}
invisible(pv)
}
pCI(Na=100,Ma=100,Mb=85)
921:卵の名無しさん
20/01/29 05:56:48.55 qmYv2xxr.net
>>868
この手法で新薬の血中濃度は老人でも統計的に有意差はないので、安心して老人にも処方できます
という商用パンプを目にする。
底辺シリツ医=裏口容疑者はこういうのに騙されるんだよなぁ。
922:卵の名無しさん
20/01/29 07:45:56.70 qmYv2xxr.net
>>868
平均値100 標準偏差15で定義される知能指数で
標準大学の新入生の知能指数の平均が100
裏口シリツ医大の新入生の知能指数の平均が85であったとする。
1学年は100人として
各大学から1/10を無作為抽出して知能指数をt検定したときのp値の期待値、中央値を求めよ。
また、p値が0.05以上になって裏口シリツ医大の新入生の知能指数は統計的に有意差はないと主張できる確率はいくらか?
923:卵の名無しさん
20/01/29 10:58:46.40 wNUvDIoi.net
>>863
大漁っ!!イサキぃぃ!!おにいちゃんかっこいいいいぃぃぃい ぃくううううう!
924:卵の名無しさん
20/01/29 12:49:44 NU71nZE3.net
私立出身医が一番多いスレであろう「開業医が集うスレ」のほうには書き込みせんのな?
925:卵の名無しさん
20/01/29 15:37:18.14 qmYv2xxr.net
>>872
書いているよ。
開業医スレでもシリツは馬鹿にされているじゃん。
スレリンク(hosp板:54番)
>>23
生データもある程度掲載されているから自分で統計処理が再現でいるな。
こういうのを考えてみると面白い。
Fisher testの結果が俺の計算と異なって気持ちが悪いままだが。
2019-nCov感染者41例中ICU収容13例(うち現喫煙者0),28例非収容(うち現喫煙者3)である。
このことから、煙草を吸うと重症化が防げると結論できるか?
926:卵の名無しさん
20/01/29 15:37:56.49 qmYv2xxr.net
>>871
そんなにレスじゃ、フェラもイサキも手に入れることはできんぞ。
927:卵の名無しさん
20/01/29 17:41:57.38 wNUvDIoi.net
では次の問題です
タバコは脳の病気です
この事実を証明せよ
URLリンク(togetter.com)
928:卵の名無しさん
20/01/29 17:46:37.77 qmYv2xxr.net
>853
929:卵の名無しさん
20/01/29 19:50:54 qmYv2xxr.net
>>872
俺が書かなくてもこういう投稿を目にする。
まあ、シリツ医という表記は俺の影響かもしれん。
自治医大とかを除外する意味で私立と書いていないわけだが。
>>
優秀な女子はぜひ国立医に入ってほしい
シリツ医は開業医のアホ子弟専用の専門学校
930:卵の名無しさん
20/01/29 21:34:38.87 wNUvDIoi.net
まあたしかに私立は学費が高いのは認める
わしも経済的理由で
合格した慶応は学費が高いので蹴って
地元イナカの国立にいった
931:卵の名無しさん
20/01/30 05:17:39 m/EdOA
932:9B.net
933:卵の名無しさん
20/01/30 08:20:09.67 b9qEKmQr.net
>>879
イナカモノにしてみれば
地元以外のところに大学が存在しているという時点で
学費以外にもアパート代や食費などの参入障壁があるんだな
934:卵の名無しさん
20/01/30 08:31:17 17lSJR5q.net
>>880
シリツ医大専門受験予備校のバイトは一コマ90分で手取り2万円だったな。大学の授業料が年間14万4千円の頃の話。
地方にはそんな割のいいバイトはない。
935:卵の名無しさん
20/01/30 20:25:51.32 m/EdOA9B.net
Rで困ったらstackoveflow.comだな。
library(Rmpfr)
vect = rep(0, 5)
for(i in 1:5){
vect[i] = mpfr(x=10^-(10*i), precBits=100)
}
# My vector has just turned to a list
vect
# Sum of list is an error
sum(vect)
# Turn it to a vector
converted_vect = Rmpfr::mpfr2array(vect, dim = length(vect))
converted_vect
# Now my sum and prod work fine and the precision is not lost
sum(converted_vect)
prod(converted_vect)
936:卵の名無しさん
20/01/30 21:10:39 m/EdOA9B.net
パッケージの内部関数を使って解決してくれた。
The function mpfr2array is not suposed to be called by the user, it is an internal tool for the package. However it's one way to solve the problem.
937:卵の名無しさん
20/01/30 21:21:54.99 m/EdOA9B.net
Rはデフォルトでは不定長さ整数を扱えないので工夫がいる。
その点はpythonやHaskellの方が有利。
"
H(n) = Σ[k=1,2,...,n] 1/k
とする。H(n)を既約分数で表したときの分子の整数をf(n)と表す。
(1)lim[n→∞] H(n) を求めよ、答えのみで良い。
(2)n=1,2,...に対して、f(n)に現れる1桁の整数を全て求めよ
"
rm(list=ls())
H <- function(n,prec=1000){ # Σ 1/kを既約分数表示する
GCD <- function(a,b){ # ユークリッドの互除法
r = a%%b # a=bq+r ⇒ a%%b=b%%rで最大公約数表示
while(r!=0){a = b ; b = r ; r = a%%b}
b }
library(Rmpfr)
one = mpfr(1, prec) # 1(one)を1000桁精度に設定
nn = 1:n # nn : 1 2 3 ... n
nume=numeric(n) # 分子の容器
for(i in nn) nume[i] = prod(nn[-i])*one # nnからi番目の要素を除いて乗算し精度アップ
nume = mpfr2array(nume, dim = length(nume)) # mpfr2arrayで加算を可能にする
Nume = sum(nume) # numeの総和を計算して分子に
Deno = factorial(n)*one # 分母のn!を精度アップ
gcd = GCD(Nume,Deno) # Numerator/Denominator約分するため最大公約数を計算
res=list(nume = Nume/gcd,deno=Deno/gcd,ratio=as.numeric(Nume/Deno)) # 最大公約数で除算して
cat(as.numeric(res$nume),'/',as.numeric(res$deno),'\n') # 数値化して分数表示
invisible(res)
}
938:卵の名無しさん
20/01/30 23:51:43 m/EdOA9B.net
#Reed-Frost Model
ReedFrost=function(
p=0.04,
N=100,
T=40)
{
q=1-p
I=numeric(T)
S=numeric(T)
I[1]=1
S[1]=N-I[1]
for(t in 1:(T-1)){
I[t+1]=S[t]*(1-q^I[t])
S[t+1]=S[t]-I[t+1]
}
plot(1:T,I,type="l",lwd=2, ylim=c(0,N
939:),xlab="time",ylab="persons",main=paste("Reed-Frost Model p= ",p)) lines(S,lty=2,col=2,lwd=2) lines(N-S,lty=3,col=3,lwd=2) legend("topright",bty="n",legend=c("Infected","Susceptible","Immunized"),lty=1:3,col=1:3,lwd=2) } par(mfrow=c(1,2)) ReedFrost(0.04) ReedFrost(0.015) ReedFrost(0.30) ReedFrost(0.03) ReedFrost(0.003)
940:卵の名無しさん
20/01/30 23:53:06 m/EdOA9B.net
こっちは、モンテカルロによるシミュレーション
# Reed-Frost and Greenwood epidemic models
# written by Dennis Chao (1/2009)
# reedfrost - the Reed-Frost epidemic model
# p = probability of transmission
# I0 = initial number of infecteds
# S0 = initial number of susceptibles
# n = number of trials
# greenwood = set to TRUE for the Greenwood model, otherwise run Reed-Frost
# outputs the number of infected and susceptibles over time (as I and S)
reedfrost <- function(p, I0, S0, n, greenwood=FALSE) {
S <- St <- rep(S0, n)
I <- It <- rep(I0, n)
q <- 1-p
time <- 0
while (sum(It)>0) {
if (greenwood)
It <- rbinom(n, St, ifelse(It>0, p, 0))
else
It <- rbinom(n, St, 1-(q^It))
St <- St-It
I <- rbind(I, It)
S <- rbind(S, St)
time <- time+1
}
I <- as.matrix(I)
S <- as.matrix(S)
list(I0=I0,S0=S0,p=p,n=n,I=I,S=S)
}
941:卵の名無しさん
20/01/31 00:59:52 24QiZJ0Y.net
f <- function(n,prec=1000){ # Σ 1/kを既約分数表示する
if(n==1){
cat(1,'\n')
invisible(1)
}else{
GCD <- function(a,b){ # ユークリッドの互除法
r = a%%b # a=bq+r ⇒ a%%b=b%%rで最大公約数表示
while(r!=0){a = b ; b = r ; r = a%%b}
b }
library(Rmpfr)
one = mpfr(1, prec) # 1(one)を1000桁精度に設定
nn = 1:n # nn : 1 2 3 ... n
nume=numeric(n) # 分子の容器
for(i in nn) nume[i] = prod(nn[-i])*one # nnからi番目の要素を除いて乗算し精度アップ
nume = mpfr2array(nume, dim = length(nume)) # mpfr2arrayで加算を可能にする
Nume = sum(nume) # numeの総和を計算して分子に
Deno=factorialZ(n) # 分母 n! = factorial(n*one)
gcd = GCD(Nume,Deno) # Numerator/Denominator約分するため最大公約数を計算
res=list(nume = Nume/gcd,deno=Deno/gcd,ratio=as.numeric(Nume/Deno)) # 最大公約数で除算して
# 分数表示 give.head=FALSEでheader除去,digits.dで桁数を指定
# capture.outputで変数に取り込み
nm = capture.output(str(res$nume, give.head=FALSE,digits.d = prec))
dn = capture.output(str(res$deno, give.head=FALSE,digits.d = prec))
cat(paste0(nm,'/',dn,'\n'))
invisible(res)
}}
for(i in 1:30) f(i)
942:卵の名無しさん
20/01/31 01:00:30 24QiZJ0Y.net
> for(i in 1:30) f(i)
1
3 / 2
11 / 6
25 / 12
137 / 60
49 / 20
363 / 140
761 / 280
7129 / 2520
7381 / 2520
83711 / 27720
86021 / 27720
1145993 / 360360
1171733 / 360360
1195757 / 360360
2436559 / 720720
42142223 / 12252240
14274301 / 4084080
275295799 / 77597520
55835135 / 15519504
18858053 / 5173168
19093197 / 5173168
122755644038509457 / 32872539188238750
186187999757029099 / 49308808782358125
14112026408124257248 / 3698160658676859375
185305423634953775872 / 48076088562799171875
5051322526706550956032 / 1298054391195577640625
11894590428248250515456 / 3028793579456347828125
1043915747995966839455744 / 263505041412702261046875
2255784105806550548873216 / 564653660170076273671875
943:卵の名無しさん
20/01/31 01:12:13.74 24QiZJ0Y.net
>>884
H(n) = Σ[k=1,2,...,n] 1/k
とする。H(n)を既約分数で表したときの分子の整数をf(n)と表す。
f(77)を求めよ。
> f(77)
17610982920730618962802441030965952272844514966520010106103127939813509744122599441432576
/ 35740194818708235
944:59764745233429885438685864430565417716720215849457565210956573486328125
945:卵の名無しさん
20/01/31 07:13:25 24QiZJ0Y.net
>>884
n=100と大きいとNAが混ざるな
946:卵の名無しさん
20/01/31 07:15:08 24QiZJ0Y.net
f <- function(n,prec=10000){ # Σ 1/kを既約分数表示する
if(n==1){
cat(? ?,1,?\n?)
invisible(1)
}else{
GCD <- function(a,b){ # ユークリッドの互除法
r = a%%b # a=bq+r ⇒ a%%b=b%%rで最大公約数表示
while(r!=0){a = b ; b = r ; r = a%%b}
b }
947:卵の名無しさん
20/01/31 07:15:14 24QiZJ0Y.net
library(Rmpfr)
one = mpfr(1, prec) # 1(one)を10000桁精度に設定
nn = 1:n # nn : 1 2 3 ... n
nume=numeric(n) # 分子の容器
for(i in nn) nume[i] = prod(nn[-i])*one # nnからi番目の要素を除いて乗算し精度アップ
nume = mpfr2array(nume, dim = length(nume)) # mpfr2arrayで加算を可能にする
Nume = sum(nume) # numeの総和を計算して分子に
Deno=factorialZ(n) # 分母 n! = factorial(n*one)
gcd = GCD(Nume,Deno) # Numerator/Denominator約分するため最大公約数を計算
res=list(nume = Nume/gcd,deno=Deno/gcd,ratio=as.numeric(Nume/Deno)) # 最大公約数で除算して
# 分数表示 give.head=FALSEでheader除去,digits.dで桁数を指定
# capture.outputで切り取りsubstrで[1]を除去
# nm = capture.output(str(res$nume, give.head=FALSE,digits.d = prec)) NAが混ざるバグあり
nm =substr(capture.output(res$nume)[2],5,nchar(res$nume)) # ?[1] 1234..? 文字列の5文字目から最後まで
dn =substr(capture.output(res$deno)[2],5,nchar(res$deno))
cat(paste0(nm,?/?,dn,?\n?))
invisible(res)
}}
948:卵の名無しさん
20/01/31 07:15:42 24QiZJ0Y.net
> f(100)
3055216077446868329553816926933899676639525195878807877583434152044192757431459126874725081455196840519615954410565802448075620352
/ 588971222367687651371627846346807888288472382883312574253249804256440585603406374176100610302040933304083276457607746124267578125
949:卵の名無しさん
20/01/31 07:45:49.87 24QiZJ0Y.net
Reed-Frost モデル
(1) 集団内の感染者と感受性のあるものとの接触はランダムに起こる
(2) 感染者と感受性のあるものが接触して伝播する確率は一定である
(3) 感染のあと必ず免疫が起こる(再感染はしない)
(4) その集団は他の集団から隔離されている
(5) 上記の条件は各時間経過中一定である
ReedFrost=function(
p=0.04, # 1期間内での伝播確率
N=100, # 集団の人数
T=40) # 全期間
{
q=1-p # 伝播させない確率
I=numeric(T) # 感染者数 Infecteds
S=numeric(T) # 感受性のある人数 Susceptibles
I[1]=1 #一人から始まったとする
S[1]=N-I[1]
for(t in 1:(T-1)){
I[t+1]=S[t]*(1-q^I[t]) # Reed-Frost Model
S[t+1]=S[t]-I[t+1]
}
return(list(I,T))
}
950:卵の名無しさん
20/01/31 11:25:54.82 24QiZJ0Y.net
# simulation model using binominal random number
rm(list=ls())
reedfrost <- function(p, I0, S0, n, greenwood=FALSE) {
S <- St <- rep(S0, n) #
951: St : Suscepibles @ time t, S: I <- It <- rep(I0, n) # It : Infected @ time t q <- 1-p # probability of non-transmission time <- 0 while (sum(It)>0) { # until no new transmission if (greenwood) It <- rbinom(n, St, ifelse(It>0, p, 0)) else It <- rbinom(n, St, 1-(q^It)) # how many ppl newly trannsmitted among St St <- St-It # how many ppl left non-infected I <- rbind(I, It) S <- rbind(S, St) time <- time+1 } Inames=NULL for(i in 0:(nrow(I)-1)) Inames[i+1]=paste0('I',i) rownames(I)=Inames Snames=NULL for(i in 0:(nrow(I)-1)) Snames[i+1]=paste0('S',i) rownames(S)=Snames list(I0=I0,S0=S0,p=p,n=n,I=I,S=S) } re=reedfrost(0.05,1,99,100) nr=nrow(re$I) plot(0:(nr-1),re$I[,1],bty='l',type='l',lwd=2, ylim=c(0,max(re$I)),xlab='time',ylab='Infecteds') for(i in 2:n) lines(0:(nr-1),re$I[,i],lwd=2,col=sample(colours()))
952:卵の名無しさん
20/01/31 15:00:17 24QiZJ0Y.net
# SEIR MODEL
"
dS(t)/dt=-bS(t)I(t),
dE(t)/dt=bS(t)I(t)-aE(t) ,
dI(t)/dt=aE(t)-gI(t) ,
dR(t)/dt=gI(t)
a:発症率,b:感染率,g:回復率
"
remove (list = objects() )
graphics.off()
SEIR <- function(
# Parameters
contact_rate = 10, # number of contacts per day
transmission_probability = 0.07, # transmission probability
infectious_period = 5, # infectious period
latent_period = 2, # latent period
# Initial values for sub-populations.
s = 9990, # susceptible hosts
e = 9, # exposed hosts
i = 1, # infectious hosts
r = 0, # recovered hosts
# Output timepoints.
timepoints = seq (0, 100, by=0.5)
){
953:卵の名無しさん
20/01/31 15:00:39 24QiZJ0Y.net
library(deSolve)
# Function to compute derivatives of the differential equations.
seir_model = function (current_timepoint, state_values, parameters)
{
# create state variables (local variables)
S = state_values [1] # susceptibles
E = state_values [2] # exposed
I = state_values [3] # infectious
R = state_values [4] # recovered
with (
as.list (parameters), # variable names within parameters can be used
{
# compute derivatives
dS = (-beta * S * I)
dE = (beta * S * I) - (delta * E)
dI = (delta * E) - (gamma * I)
dR = (gamma * I)
# combine results
results = c (dS, dE, dI, dR)
list (results)
}
)
}
# Compute values of beta (tranmission rate) and gamma (recovery rate).
beta_value = contact_rate * transmission_probability
gamma_value = 1 / infectious_period
delta_value = 1 / latent_period
954:卵の名無しさん
20/01/31 15:00:45 24QiZJ0Y.net
# Compute Ro - Reproductive number.
Ro = beta_value / gamma_value
# Disease dynamics parameters.
parameter_list = c (beta = beta_value, gamma = gamma_value, delta = delta_value)
# Compute total population.
N = s + i + r + e
# Initial state values for the differential equations.
initial_values = c (S = s/N, E = e/N, I = i/N, R = r/N)
# Simulate the SEIR epidemic.
# ?lsoda # Solver for Ordinary Differential Equations (ODE), Switching Automatically Between Stiff and Non-stiff Methods
output = lsoda (initial_values, timepoints, seir_model, parameter_list)
# Plot dynamics of Susceptibles sub-population.
plot (S ~ time, data = output, pch='S', type='b', col = 'blue', bty='l', cex=0.75,
ylab='SEIR')
points(E ~ time, data = output, pch='E', type='b', col = 'yellow', cex=0.75)
points(I ~ time, data = output, pch='I', type='b', col = 'red', cex=0.75)
points(R ~ time, data = output, pch='R', type='b', col = 'green', cex=0.75)
}
SEIR()
args(SEIR)
SEIR(contact_rate=0.5,transmission_probability=0.05,infectious_period=365,latent_period=7,s=99,e=0,i=1,r=0,timepoints = 0:365*6)
955:卵の名無しさん
20/01/31 17:42:44 24QiZJ0Y.net
エンデミックな定常状態を(S?, I?)とおけば、S?N=1R0,I?N=μμ+γ(1?1R0)(15)である。すなわちエンデミックな状態における感受性人口比率と基本再生産数は逆数関係にあり、有病率(prevalence)I?/Nは1?1/R0に比例していて、その比例係数は、感染状態における平均滞在時間1/(μ+γ)とホストの寿命1/μの比である。これらの式はエンデミックな感染症におけるR0の推定式ともみなせる
956:卵の名無しさん
20/01/31 17:45:07.01 24QiZJ0Y.net
Rv?1という条件はワクチン接種率の条件として書き直せばv?1?1R0=H
957:卵の名無しさん
20/02/01 15:12:17 hIisy8jC.net
H(n) = Σ[k=1,2,...,n] 1/k
とする。H(n)を既約分数で表したときの分子の整数をf(n)と表す。
(1)lim[n→∞] H(n) を求めよ、答えのみで良い。
(2)n=1,2,...に対して、f(n)に現れる1桁の整数を全て求めよ
"
H <- function(n) sum(1/(1:n))
plot(sapply(1:1000,H),bty='l')
f <- function(n,prec=10000){ # Σ 1/kを既約分数表示する n>>=23で誤計算
if(n==1){
cat(n, ':' ,1,'\n')
invisible(1)
}else{
GCD <- function(a,b){ # ユークリッドの互除法
r = a%%b # a=bq+r ⇒ a%%b=b%%rで最大公約数表示
while(r!=0){a = b ; b = r ; r = a%%b}
b }
nn = 1:n # nn : 1 2 3 ... n
nume=list() # 分子の容器
length(nume)=n
Nume=0
library(gmp)
for(i in nn) Nume = Nume + prod.bigz(nn[-i]) # nnからi番目の要素を除いて乗算して総和を分子に
Deno=factorialZ(n) # 分母
gcd = GCD(Nume,Deno) # Numerator/Denominator約分するため最大公約数を計算
NUME=Nume/gcd
DENO=Deno/gcd
ratio=as.numeric(Nume/Deno)
RE=list(NUME,DENO,ratio) # 最大公約数で除算して約分
cat(n, ':' ,as.character(NUME),'/',as.character(DENO),'\n')
invisible(RE)
}}
958:卵の名無しさん
20/02/01 22:05:31.53 hIisy8jC.net
H(n) = Σ[k=1,2,...,n] 1/k
とする。H(n)を既約分数で表したときの分子の整数をf(n)と表す。
f <- function(n){ # Σ 1/kを既約分数表示する
if(n==1){
cat(n, ':' ,1,'\n')
invisible(1)
}else{
GCD <- function(a,b){ # ユークリッドの互除法
r = a%%b # a=bq+r ⇒ a%%b=b%%rで最大公約数表示
while(r!=0){a = b ; b = r ; r = a%%b}
b }
nn = 1:n # nn : 1 2 3 ... n
nume=list() # 分子の容器となるlist
length(nume)=n # を設定
Nume=0
library(gmp)
for(i in nn) Nume = Nume + prod.bigz(nn[-i]) # nnからi番目の要素を除いて乗算して総和を分子に
Deno=factorialZ(n) # 分母
gcd = GCD(Nume,Deno) # 約分するため最大公約数を計算
NUME=Nume/gcd
DENO=Deno/gcd
ratio=as.numeric(Nume/Deno)
RE=list(NUME,DENO,ratio)
cat(n, ':' ,as.character(NUME),'/',as.character(DENO),'\n')
invisible(RE)
}}
959:卵の名無しさん
20/02/02 13:53:23 +5dNqMpE.net
tbl <- function(x,v){ # vの要素がxにいくつあるか集計する
n=length(v)
hme=numeric(n) # how many entries?
for(i in 1:n) hme[i]=sum(x==v[i])
rbind(v,hme)
}
tbl(sample(10,rep=T),1:10)
960:卵の名無しさん
20/02/04 16:35:57 b64IHQrg.net
"URLリンク(this.kiji.is)
中国湖北省武漢市からチャーター機で日本へ帰国した邦人の新型コロナウイルス感染率が高いと、
中国で驚きの声が上がっている。中国当局が発表した同市の感染者の割合に比べ「39倍も高い」というのだ。
現地は医療現場が混乱しているため、実際には発表よりかなり多くの感染者がいる可能性がある。
日本政府はチャーター機計3便を武漢市に派遣し、邦人565人が帰国した。厚生労働省によると、
チャーター機に乗っていた感染者は、症状のない人も含め計8人。感染率は1.416%だ。
一方、1月31日現在、武漢市の感染者数は3215人で、感染率は0.036%にとどまった。"
r1=8
r2=3215
n1=565
n2=round(3215/(0.036/100))
poisson.test(c(r1,r2),c(n1,n2))
prop.test(c(r1,r2),c(n1,n2))
library(BayesFactor)
mat=cbind(c(r1,r2),c(n1,n2)-c(r1,r2))
fisher.test(mat)
chisq.test(mat)$p.value
bf=BayesFactor::contingencyTableBF(mat,sampleType = 'poisson',posterior = TRUE,iteration=1e5)
rbf=round(bf)
(r1/n1)/(r2/n2)
frp <- function(x) (x[1]/(x[1]+x[2]))/ (x[3]/(x[3]+x[4]))
rrp=apply(rbf,1,frp)
hist(rrp)
BEST::plotPost(rrp)
(r1/(n1-r1))/(r2/(n2-r2))
fop <- function(x) (x[1]/x[2])/ (x[3]/x[4])
orp=apply(rbf,1,fop)
BEST::plotPost(orp)
quantile(orp,c(0.005,0.025,0.5,0.975,0.995))
library(ex
961:act2x2) exact2x2::fisher.exact(mat) fisher.test(mat)
962:卵の名無しさん
20/02/04 16:57:56 b64IHQrg.net
# Reed-Frost Model 時刻t+1での感染者数=時刻tでの感受性人数*(少なくとも一人の感染者がでる確率)
# p=0.04;N=100;I0=1;T=8
ReedFrost=function(
p=8/565, # 1期間内での伝播確率
N=565, # 集団の人数
I0=8, # 最初の感染者数
T=30, # 全期間
print=TRUE # グラフ表示
)
{
q=1-p # 伝播させない確率
I=numeric(T) # 感染者数 newly Infecteds
S=numeric(T) # 感受性のある人数 Susceptibles
I[1]=I0 # I0人から始まったとする
S[1]=N-I[1]
for(t in 1:(T-1)){
I[t+1]=S[t]*(1-q^I[t]) # Reed-Frost Model
S[t+1]=S[t]-I[t+1]
}
if(print){
plot(1:T,I,type='b',pch=19,lwd=2, ylim=c(0,N),bty='l',
xlab="time",ylab="persons",main=paste("Reed-Frost Model p= ",round(p,3)))
points(S,lty=2,col=2,lwd=2,type='b',pch=19)
points(N-S,lty=3,col=3,lwd=2,type='b',pch=19)
legend("right",bty="n",legend=c("New Infection","Susceptible","Immunized"),lty=1:3,col=1:3,lwd=2)}
list(newly_infecteds=round(I),susceptibles=round(S))
}
ReedFrost(T=7)
963:卵の名無しさん
20/02/06 22:51:49 WriFO9uC.net
"厚生労働省によりますと、新型コロナウイルスに感染した香港の男性が乗っていたクルーズ船で、
症状がみられる人やその濃厚接触者などにウイルス検査を実施した結果、10人が感染していたことが新たに確認されました。
これまでにこのクルーズ船の乗員乗客で感染が確認されたのは、合わせて20人です。
厚生労働省によりますと、クルーズ船「ダイヤモンド・プリンセス」の船内では
乗客と乗組員全員の合わせて3700人余りの検疫が行われ、発熱やせきなどの症状があった120人と、
症状がある人と濃厚接触した153人の合わせて273人から検体を採取して順次、ウイルス検査を実施しています。
6日、新たに71人分の結果が判明しこのうち10人がウイルスに感染していたことが確認されました。
"
NN=3700
N=273
n=71
r=10
(re=binom::binom.confint(r,n,conf=0.95))
cbind(re[,1],re[,5:6]*N)
library(BayesFactor)
(mat=cbind(c(r,r),c(n-r,n-r)))
bf=BayesFactor::proportionBF(r,n,r/n,posterior = TRUE,iterations=10000)
plot(bf)
head(bf)
infct=bf[,'p']*N
summary(infct)
BEST::plotPost(N*bf[,2],credMass = 0.95)
964:卵の名無しさん
20/02/12 15:09:27 Th6rMqDv.net
んほおぉぉ
URLリンク(m.youtube.com)
965:卵の名無しさん
20/02/12 15:41:37.89 6n8f5vvH.net
日本最高学費の底辺私立医大では
1年:進級失敗10人
2年:進級失敗16人
3年:進級失敗34人
4年:進級失敗9人
5年:進級失敗10人
6年:卒業失敗26人
一学年約120~130人前後。
同じ学年で二回留年すると退学
スレリンク(doctor板:1番)
であるという。
1年次学費総額 12,145,000円 2年次以降学費(年間) 7,030,000円
1学年を125人として上記データから算出した確率(例、1年次は10/125の確率で留年)を用いて
卒業できる確率と卒業生の在学年数の期待値を求めよ。
また、退学になる確率と退学者の在学年数の期待値を求めよ。
966:卵の名無しさん
20/02/12 15:41:57.91 6n8f5vvH.net
(p=1-c(10,16,34,9,10,26)/125)
q=1-p # 留年確率
(P=prod(1-q^2)) # 卒業できる確率 Π( 1 - 2年連続留年確率)
(Q=1-P) # 退学となる確率
p1=p # 各学年を1年で進級する確率
p2=p*(1-p) # 各学年を2年で進級する確率
P12=rbind(p1,p2)
gr=expand.grid(1:2,1:2,1:2,1:2,1:2,1:2) # 6学年の組み合わせdata.frame
years=apply(gr,1,sum) # 各組合せでの在学年数
dat=cbind(gr,years)
dat=as.matrix(dat) # data.frameを行列化
f <- function(x){ # 1~6年次の進級に
967:かかった年数の配列xからその確率を返す ps=1 # 進学確率 for(i in 1:6){ # 各組合せでの通算確率を算出 ps=ps*P12[x[i],i] } return(ps) # 卒業確率(入学してからその組合せで卒業できる確率) } Ps=apply(dat,1,f) # 組合わせ候補にfを適用 Probability of success data=cbind(dat,Ps) # dat[,7]:入学してからの確率 dat[,8]:在学年数 # 在学年数の期待値 sum(data[,7]*data[,8])/sum(Ps) # 期待値=Σxi*PiでΣPi=1でなくてはならないのでsum(Ps)で除算 # 退学 q=c(10,16,34,9,10,26)/125 Q=numeric(6) Q[1]=q[1]^2 # 1年次で退学 Q[2]=q[2]^2 *(1- Q[1]) # 2年次で退学 for(i in 2:6){ Q[i]=q[i]^2 * (1-Q[i-1]) } sum((1:6)*Q)/sum(Q)
968:卵の名無しさん
20/02/12 15:42:45.68 6n8f5vvH.net
ko=vector('list',6) # kicked out 退学までの軌跡
a=1:2 # 1:1回で進級,2:1回留年して翌年進級
ko[[1]]=as.matrix(2) # 1年時で退学
ko[[2]]=as.matrix(expand.grid(a,2)) # 2年時で退学の軌跡
ko[[3]]=as.matrix(expand.grid(a,a,2)) # 3年次で退学の軌跡
ko[[4]]=as.matrix(expand.grid(a,a,a,2))
ko[[5]]=as.matrix(expand.grid(a,a,a,a,2))
ko[[6]]=as.matrix(expand.grid(a,a,a,a,a,2))
f <- function(x){ # 軌跡与えてその確率(入学からの確率なので総和が退学確率Qに等しい)
n=length(x)
if(n==1) return(q[1]^2)
Pf=1
for(i in 1:(n-1)) Pf=Pf*P12[x[i],i]
return(Pf*q[n]^2)
}
Ef=numeric(6) # i年次退学する時の在学年数期待値(Qで除算が必要)
for(i in 1:6){
pf=apply(ko[[i]],1,f) # i年次で退学する場合の軌跡での確率
nf=apply(ko[[i]],1,sum) # その軌跡での在学年数
Ef[i]=sum(pf*nf)
}
sum(Ef)/Q
969:卵の名無しさん
20/02/12 15:43:16.39 6n8f5vvH.net
日本最高学費の底辺私立医大では
1年:進級失敗10人
2年:進級失敗16人
3年:進級失敗34人
4年:進級失敗9人
5年:進級失敗10人
6年:卒業失敗26人
一学年約120~130人前後。
同じ学年で二回留年すると退学
スレリンク(doctor板:1番)
であるという。
1年次学費総額 12,145,000円 2年次以降学費(年間) 7,030,000円
父兄から 同じ学年で二回留年すると退学というルールは厳しすぎる、
これを撤廃して、最長12年までは在学可能とすべきだという提言がなされたとする
この提言を採用した場合、
1学年を125人として上記データから算出した確率(例、1年次は10/125の確率で留年)を用いて
卒業できる確率と卒業生の在学年数の期待値を求めよ。
また、退学になる確率と退学者の在学年数の期待値を求めよ。
970:、
20/02/12 20:50:03 dWksx8Sa.net
こんな凄いスレッドあったんですね、読み返してみます、
971:卵の名無しさん
20/02/12 21:02:12.30 iCMxhp27.net
ここは国試浪人の事務員のスレ
972:卵の名無しさん
20/02/12 21:21:16 6n8f5vvH.net
>>913
手コキ問題の答はまだか?
973:卵の名無しさん
20/02/12 22:17:18 6n8f5vvH.net
>>911
まず、シミュレーションプログラムを作る。
p=1-c(10,16,34,9,10,26)/125 #進級確率
q=1-p
sim <- function(){
g=1 # 学年,7:卒業
y=0 # 在学年数
while(g<7 & y<12){
g=g+sample(1:0,1,prob=c(p[g],q[g]))
y=y+1
}
c(g,y)
}
k=1e3
replicate(k,,sim())
974:卵の名無しさん
20/02/13 13:57:25 DpKqMtsQ.net
p=1-c(10,16,34,9,10,26)/125 #進級確率
q=1-p
sim <- function(){
g=1 # 学年,7:卒業
y=0 # 在学年数
while(g<7 & y<12){ # 卒業前、在学12年未満なら
g=g+sample(1:0,1,prob=c(p[g],q[g])) # 進級判定して
y=y+1 # 在学年数を増やす
}
c(g,y) # 学年と在学年数を返す
}
k=1e6
RE=t(replicate(k,sim()))
head(RE)
summary(RE)
R12=RE[RE[,2]==12,] # 12年在学
nrow(R12[R12[,1]<7,])/k # 退学確率
975:卵の名無しさん
20/02/13 13:58:30 DpKqMtsQ.net
> nrow(R12[R12[,1]<7,])/k # 退学確率
[1] 0.001039
976:卵の名無しさん
20/02/13 21:06:00 DpKqMtsQ.net
p=1-c(10,16,34,9,10,26)/125 # 進級確率
q=1-p # 留年確率
p6=prod(p) # 6年で卒業確率
f <- function(x) prod(q[x])
pf=c(6,p6)
for(i in 7:12){
cmb=combinations(6,i-6,rep=T)
pf=rbind(pf,c(i,p6*sum(apply(cmb,1,f))))
}
pf
(P=sum(pf[,2]))# 12年以内に卒業できる確率
sum(pf[,1]*pf[,2])/P # 卒業までの平均在学年数
977:卵の名無しさん
20/02/14 01:05:55.44 9rZr8ZNE.net
library(gtools)
p=1-c(10,16,34,9,10,26)/125 # 進級確率
q=1-p # 留年確率
p6=prod(p) # 6年で卒業確率
f <- function(x) prod(q[x]) # x:留年学年の配列→その組合せでの確率
pf=c(6,p6) # 在学年数,その確率
for(i in 7:12){#在学年数i 7~12において
cmb=combinations(6,i-6,rep=T) # 1:6から留年学年の組合せ
978:i-6個を重複を許して選ぶ pf=rbind(pf,c(i,p6*sum(apply(cmb,1,f)))) #その組合せでの確率の総和に6年分の進級確率を乗ずる } (P=sum(pf[,2]))# 12年以内に卒業できる確率 (yg=sum(pf[,1]*pf[,2])/P) # 卒業までの平均在学年数 Σ x*pdf (Q=1-P) # 退学確率 P*yg+Q*12 # # 大学を出る(卒業もしくは退学)を平均年数
979:卵の名無しさん
20/02/14 01:13:37.89 9rZr8ZNE.net
理論値
> (yg=sum(pf[,1]*pf[,2])/P) # 卒業までの平均在学年数 Σ x*pdf
[1] 7.027985
シミュレーション値
> mean(R7[,2]) # 卒業生の在学年数
[1] 7.027643
合致しているから、どちらも正しそう。
980:卵の名無しさん
20/02/14 18:43:07 x9TXrlfR.net
思えば東京に住んでる奴はそれだけで経済的に優遇されてるな
田舎に住んでる貧乏人は学力があっても地元国立一択だったり
学力があっても就職コースだったりするからなぁ
981:卵の名無しさん
20/02/14 20:30:31.44 9rZr8ZNE.net
>>921
東京で地元の国立に進学するのは容易ではないと思う。
982:卵の名無しさん
20/02/14 20:33:12.20 9rZr8ZNE.net
>>919
この問題は数学板に投稿してもプログラムの練習にしかならんな。
983:卵の名無しさん
20/02/14 22:03:29.56 wJciaHaO.net
難し、ネットで医療統計のページ見つけたが、
エクセルシートが出てくるだけだし、χ2検定だけでも
俺みたいなのに判るように解説してほしい、
MRが言う当社の薬剤は他のよりも優れていない確率は5%
って言うのはホントなんだろうか?違うよね
984:卵の名無しさん
20/02/15 07:37:08.29 6ICubIld.net
>>924
Χ2乗検定が知りたいのか
p値で判定することの欺瞞性を知りたいのか
非劣性検定の問題点を知りたいのか
どれ?
985:卵の名無しさん
20/02/15 15:24:15 BqM3deVD.net
具体的にはグーフィス(持田)の製品紹介にある
グーフィスはプラセボ群に比して有意に大きな値をしめしました(p<0.0001共分散解析)
次行にはほぼ同じ文面で同じ数字でFisherの正確検定
更にしたには Wilconxonの順位和検定とあります
なんのことやら、
986:卵の名無しさん
20/02/15 16:05:12.72 9KqLKmld.net
分散分析は郡内分散と群間分散の比がF分布に従うことをつかった検定。
(この比をF ratioというのが色っぽいw)
統計の教科書読めばFisherもWilcoxonも載っているから独習してください。
頻度主義統計を理解する前に、
p値で判定することの欺瞞性や非劣性検定の問題点を論じても混乱するだけだと思います。
987:卵の名無しさん
20/02/15 19:46:00 9KqLKmld.net
新型コロナウイルスが新たに67人で陽性であった。
検査可能数を100~1000と仮定したときの検査数の期待値と95%信頼区間を求めよ
"
f <- function(N,conf=0.95,print=FALSE){
# N=100 ; conf=0.95
m=1
r=67
n=r:N
pmf=choose(r-1,m-1)/choose(n,m) #Pr(max=60|n)
# plot(n,pmf,ylab='probability',bty='l')
pdf=pmf/sum (pmf)
if(print) plot(n,pdf,ylab='density',type='h',bty='l')
E=sum(n*pdf) #E(n)
cdf=cumsum(pdf)
Lidx=which.min(abs(cdf-(1-conf)/2))
Uidx=which.min(abs(cdf-(1-(1-conf)/2)))
return(c(E=E,L=n[Lidx],U=n[Uidx]))
}
N=seq(100,1000,by=50)
RE=as.matrix(sapply(N,f))
RE
plot(N,RE[1,],bty='l',pch=19,ylim=c(0,1000),xlab='上限件数',ylab='推定件数')
segments(x0=N,y0=RE[2,],x1=N,y1=RE[3,])
988:卵の名無しさん
20/02/16 08:15:55 dI8Llzzn.net
コロナ患者はおことわり
当たり前だよなぁ
普通の病院に来られても検査はできん
医師なら当然知ってることだ
989:卵の名無しさん
20/02/18 12:09:44 OUmetWs8.net
製薬会社の社員って基本的に自社製品の宣伝しかしないよね
金儲けのため病院訪問してんだから当たり前かも知らんが
コンプ薬屋クンとかサラリーマン時代はどんな営業してたのかねぇ?
990:卵の名無しさん
20/02/18 19:48:55 OUmetWs8.net
人生に失敗してる奴ほど反日野党に肩入れするよねW
991:卵の名無しさん
20/02/18 19:53:12 Zlyy45UZ.net
binom::binom.bayes(4,4,prior.shape1=1,prior.shape2 = 1)
curve(dbeta(x,5,1),bty='l')
options(digits=22)
pbeta(0.01,5,1)
binom::binom.bayes(4,4)
curve(dbeta(x,4.5,0.5),bty='l')
pbeta(0.01,4.5,0.5)
992:卵の名無しさん
20/02/18 19:54:51 Zlyy45UZ.net
トップ為政者がどんなに私的な放蕩や汚職等の悪行を重ねても、国の危機には民の命を第一に守るという最低限の義務を放棄しない限りは支持され続ける、アベはそれをあっさりと放棄した。
993:卵の名無しさん
20/02/18 20:07:46 Zlyy45UZ.net
日本の領土、日本民族、そして日本語をもっとも破壊している反日勢力が安倍なんだが。
994:卵の名無しさん
20/02/19 13:39:37.63 PmndJ4fY.net
はいはい ワロスワロスW
995:卵の名無しさん
20/02/19 13:51:46 JGE03rN8.net
危機に対する面構えのこの違い
URLリンク(Imgur.com)
URLリンク(Imgur.com)
これだけで一目瞭然だろう
996:卵の名無しさん
20/02/19 13:52:46 JGE03rN8.net
【コロナ】安倍内閣がクソすぎて無関心層が目覚め出した…「安倍総理ってもしかして無能なの…?」★19
スレリンク(newsplus板)
997:卵の名無しさん
20/02/20 18:44:21 i8qqM3yt.net
お兄ちゃん、イサキは釣れたのっ?
998:卵の名無しさん
20/02/24 08:53:05.25 TQ/flAIX.net
おーい国試浪人 および コンプ薬屋
お前らのスレはここだ
他スレを荒らすんじゃないぞ
999:卵の名無しさん
20/02/24 19:15:49.10 x/xpTJWO.net
一辺10[m]の正方形ABCDのプールがある
点Dでは水が湧き出しており、点Dからr[m]離れた場所では(r/10)[m/s]までのスピードでしか泳げない
点Aから正方形の中心まで泳ぐのに掛かる最短時間を求めよ
極座標の曲線r=f(θ), 10=f(-π/2), 5√2=f(-π/4)上を泳ぐ時間は
T(f)=10∫[-π/2,-π/4]√(1+(f'(θ)/f(θ))^2)dθ
δT(f)=0に対するオイラーラグランジュの方程式は
-(f'/f)^2/(f√(1+(f'/f)^2))-(d/dθ)((f'/f)/(f√(1+(f'/f)^2)))=0
整理すると
(f'^2-f''f)(f^2+f'^2)^(-3/2)=0
この解はf(θ)=a e^(bθ)で境界条件を合わせると
f(θ)=10e^((-θ-π/2)(2log2)/π)
このとき
T(f)=10∫[-π/2,-π/4]√(1+(2log2/π)^2)/dθ
=(5/2)√(π^2+4(log2)^2)
1000:卵の名無しさん
20/02/24 19:54:31.51 qbJ7+ra+.net
では、それを、みさくら語でどうぞ
1001:卵の名無しさん
20/02/24 21
1002::31:49.62 ID:x/xpTJWO.net
1003:卵の名無しさん
20/02/24 21:57:22.22 iv/jnKeV.net
底辺の人間ほどエセの正義を語ることで
偉くなった気分になるんだよね
URLリンク(togetter.com)
1004:卵の名無しさん
20/02/24 22:06:34.38 x/xpTJWO.net
Aチーム6人 : Bチーム6人でレースゲームを16試合する。
1位15点
2位12点
3位10点
4位9点
5位8点
6位7点
7位6点
8位5点
9位4点
10位3点
11位2点
12位1点
の配点の場合、各チームの合計点が656:656の引き分けになる確率を教えてください。
x=c(1:10,12,15)
sim <- function() sum(replicate(16, sum(x[sample(12,6)])))==656
k=1e5
mean(replicate(k,sim()))
1005:卵の名無しさん
20/02/24 22:18:39.05 x/xpTJWO.net
3861707060011302197274473352442662107544339496/924^16
=0.0136781633711920174806098975...
シミュレーションでの近似値が
> mean(replicate(k,sim()))
[1] 0.01344
大体あってる。
1006:卵の名無しさん
20/02/24 22:26:32.20 x/xpTJWO.net
>>941
数くらい数えるのと、足し算くらいできるだろ?
Aチーム6人 : Bチーム6人でレースゲームを16試合する。
1位15点
2位12点
3位10点
4位9点
5位8点
6位7点
7位6点
8位5点
9位4点
10位3点
11位2点
12位1点
の配点の場合、各チームの合計点が656:656の引き分けになる
282326431934901209944756971232412938108921708544通りの中から引き分けになるのが何個あるか数えてくれ。
1007:卵の名無しさん
20/02/25 22:16:10.30 dZGWYWaM.net
臨床経験が無いのに
臨床統計を語るのは
どういうお気持ちですか?
1008:卵の名無しさん
20/02/26 06:27:56.22 XNNbo5EQ.net
ここをコンプ薬屋と事務員の政治スレにしよう
1009:卵の名無しさん
20/02/26 06:41:37.50 P+li/ooa.net
>>947
臨床でのエピソード
23 卵の名無しさん sage 2020/02/06(木) 10:18:40.39 ID:fvuyjOo1
予約の患者が朝食とったらしくキャンセルで空き時間ができた。
余談ながら、
若い看護師が若い女性に(キイロカインを)「ごっくんしてください」というのは艶かしいな。
思わずVater乳頭の観察に時間を割いてしまうw
1010:卵の名無しさん
20/02/26 07:03:20.27 P+li/ooa.net
>>948
安倍ちゃんのいいところを書けば250円貰えるぞwww
://i.imgur.com/jfBeNNK.jpg
1011:卵の名無しさん
20/02/26 13:42:04.22 AcA2+/yJ.net
遥か昔の学生時代の記憶はどうでもいいから
早くコンプ薬屋も呼んできて
二人で壁打ちしてろ
1012:卵の名無しさん
20/02/26 14:21:00.17 P+li/ooa.net
>>951
また、逃げるの?
臨床医学って確率的事象を扱う稼業なのになぁ。
自分で統計処理できないのは恥ずかしいなぁ。
前スレでも答えられずに逃亡していた。
【新型コロナ】中国の映画監督、新型コロナウイルスで一家4人死亡
スレリンク(newsplus板)
上記をみた統計できない患者から、かかったら必ず死ぬんですか?
と�
1013:キかれたら数字を出して説明したいよね。 4人感染して4人死亡。 死亡率の期待値とその95%信頼区間は?
1014:卵の名無しさん
20/02/26 14:38:28.03 P+li/ooa.net
"
医学雑誌ランセットには、武漢市金銀潭医院に入院していた重度の肺炎患者52人のうち32人(61.5%)が死亡したという論文が発表されています。
URLリンク(news.yahoo.co.jp)
厚生労働省によりますと、新型コロナウイルスへの感染が確認された人で人工呼吸器をつけたり集中治療室で治療を受けたりしている重症者は、25日の時点でクルーズ船の乗客・乗員が36人、国内で感染が確認された人が14人で合わせて50人となっています。
URLリンク(www3.nhk.or.jp)
何をもって重症としているのかは判然としないが、両者の重症判定基準が同等として
50人中何人死亡すると推定されるかを95%信頼区間とともに述べよ。
"
library(binom)
(bn=binom.confint(32,52))
attach(bn)
data.frame(method,mean*50,lower*50,upper*50)
1015:卵の名無しさん
20/02/26 16:12:56.43 E0xdg3nZ.net
>>953
童貞! 校庭! 童貞! 校庭!
童貞! 校庭! 童貞! 校庭!
1016:卵の名無しさん
20/02/27 18:46:44 N3R7wpqC.net
>>953
臨床童貞はおもしろいですか?
1017:卵の名無しさん
20/02/28 12:33:39.16 Fl09KmOa.net
ド底辺シリツ医および同等頭脳が答が出せない問題。
若い女性研修医(嘘つきかどうかは不明)から
「あなたのいうことが正しければ手コキかフェラをしてあげる、そうでなければセンズリを命じる」と言われた。
フェラをしてもらうには何と言えばいいか?
尚、正直者なら嘘をつかないし、嘘つきなら必ず嘘をつく。
1018:卵の名無しさん
20/02/28 18:59:27.82 SCfPKsWW.net
>>956
いつもどおりキモいですねW
1019:卵の名無しさん
20/02/28 19:15:37 04W3JdR8.net
>>957
ド底辺シリツ医頭脳登場w
1020:卵の名無しさん
20/02/28 21:16:39 x0SPnMpa.net
>>958
早くド底辺仲間のコンプ薬屋を呼んできなよ
オマエとコンプが揃ってAGWPする
地獄絵図が見たいなぁ
by 都内JK
1021:卵の名無しさん
20/02/29 03:33:39.80 Ix6SlPJQ.net
>>959
あんたが>956に答が出せないド底辺頭脳ではないかよ。
1022:卵の名無しさん
20/02/29 07:34:46 1jppDb+T.net
>>960
コンプ薬屋とAWPまだぁ
by 都内JK
1023:卵の名無しさん
20/02/29 07:39:32.62 Ix6SlPJQ.net
>>959
>956の答も出せない底辺頭脳を晒して嬉しい?
1024:卵の名無しさん
20/02/29 10:28:46 Ix6SlPJQ.net
プログラム解はこのスレに既出だからそれを日常言語化すれば答になるけど、ド底辺はそういう教養を欠いているんだよなぁ。
1025:卵の名無しさん
20/02/29 17:53:19.62 1jppDb+T.net
相変わらず、反日テレビ局は偏向報道が酷いな
まるで何もかも反対しかしない
昭和の野党みたいだ
1026:卵の名無しさん
20/03/01 15:18:41 tomslBB/.net
おい統計野郎
他スレでのオマエのゴミレスみると
ごまかしと暇つぶしだけの人生だなW
1027:卵の名無しさん
20/03/01 15:27:55 8WrH5vXE.net
レベルが高いからレスがつかんのよ
あんたこういうグラフ作れる?
暇つぶしに、感度80% 特異度90%だった時の有病率と 検査陰性での感染確率のグラフを書いてみた。
1028: https://i.imgur.com/yVwE5Bs.png
1029:卵の名無しさん
20/03/01 15:45:45 8WrH5vXE.net
>>965
数学板にも投稿してますが、みてないの?
357 132人目の素数さん sage 2020/02/29(土) 20:50:00.58 ID:duevA7i1
>>345
RでSEIR MODEL
dS(t)/dt = mu*(N-S) - b*S(t)*I(t)/N - nu*S(t)
dE(t)/dt = b*S(t)I(t)/N - (mu+sig)*E(t)
dI(t)/dt = sig*E(t) - (mu+g)*I(t)
dR(t)/dt = g*I(t) - mu*R + nu*S(t)
mu:自然死亡率 b:感染率(S->I)
nu:ワクチン有効率(S->R) sig:発症率(E->I),g:回復率(I->R)
でプログラムを組んでクルーズ船に閉じ込めておいたときの収束予想をだそうと遊んでみた。
1030:卵の名無しさん
20/03/02 21:47:53.33 /IzIBX47.net
ほらまたコンプがあっちで湧いた
お前の嫁だろ
早く呼んで来いやー
1031:卵の名無しさん
20/03/02 22:05:34 BcEADux3.net
>>968
>956の答も出せないドアホ登場。
ド底辺シリツ医並の頭脳だなぁ。
1032:卵の名無しさん
20/03/02 23:34:26.03 /IzIBX47.net
>>969
このFランって オマエのことだろ
URLリンク(www.nicovideo.jp)
私はあなたと違って高学歴ですが
1033:卵の名無しさん
20/03/02 23:53:44.04 JuR2B+30.net
息を吐くようにシリツの言葉が出るって、どう考えても底辺国立だな
1034:卵の名無しさん
20/03/03 06:12:16.13 MOd/ilOT.net
>>971
都内旧二期校卒
あんたは>956の答も出せないドアホだろ
ド底辺シリツ医並の頭脳だなぁ。
1035:卵の名無しさん
20/03/03 06:20:52.12 MOd/ilOT.net
ド底辺シリツ医頭脳じゃなければ>956に正答すると思っていた。
それを受けてこっちをやってみ、と書く予定だったが
想定したよりも馬鹿だった。
AからHの8人はそれぞれ正直者か嘘つきであり、誰が正直者か嘘つきかはお互いに知っている。
A,B,C,D,Eは嘘つきなら必ず嘘をつくが、F,G,Hは嘘つきでも正しいことを言う場合がある。
次の証言から誰を確実に正直者と断定できるか?
A「嘘つきの方が正直者より多い」
B「Hは嘘つきである」
C「Bは嘘つきである」
D「CもFも嘘つきである」
E「8人の中に、少なくとも1人嘘つきがいる」
F「8人の中に、少なくとも2人嘘つきがいる」
G「Eは嘘つきである」
H「AもFも正直者である」
1036:卵の名無しさん
20/03/03 07:06:47.46 9I6dVjmr.net
>>972
相変わらず毎日ウソばっかりですね
このFランの恐ろしいところは周囲にうつそうとすることですね
高学歴の私は抗体があるので無事ですが
URLリンク(www.nicovideo.jp)
1037:卵の名無しさん
20/03/03 07:49:46.32 MOd/ilOT.net
>>974
あんたは>956の答も出せないドアホだろ
ド底辺シリツ医並の頭脳だなぁ。
1038:卵の名無しさん
20/03/03 07:59:50.55 MOd/ilOT.net
>>974
>973に論理的に答られたら高学歴が実証されるからやってみ。
24時間待ってあげる。
1039:卵の名無しさん
20/03/03 10:05:44.33 paZpS0zc.net
都内旧二期校卒!?とか言う言葉が出るあたり、相当な負け犬根性が培われてるな
一期校より下、旧六より下、都内なら私立より下まではいかんが私立より上ではない
底辺国立でいいだろ、分かりやすくて
高校の時に頑張って勉強しとくんだったなwww
ま、そんなことより、リアルでもネットでも、他人とコミュニケーション取れるようにな
1040:卵の名無しさん
20/03/03 16:24:05 MOd/ilOT.net
>>977
>973できた?
1041:卵の名無しさん
20/03/04 17:06:06.71 0bS4J8i5.net
>>978
受験を控えている高校生です
他のスレであなたの通っていた学校は実はFランであるとお聞きしました
レスも私の知っている大人の方たちと比べて随分と支離滅裂で低レベルなものしかされてませんが
今はどのようなお仕事をされてるんですか?
1042:卵の名無しさん
20/03/04 17:23:24.90 gP3NYsfY.net
>>979
スレリンク(hosp板:25番)
1043:卵の名無しさん
20/03/04 17:24:09.67 gP3NYsfY.net
>>979
>973に代わり答えてみ!
1044:卵の名無しさん
20/03/04 17:57:44.47 JeuwPZj2.net
>>981
それで今はどのようなお仕事をされてるんですか?
下品なアジを書いたプラカード持って街角に立ってるお仕事ですか?
1045:卵の名無しさん
20/03/04 18:02:57.24 gP3NYsfY.net
>>982
スレリンク(hosp板:25番)
1046:卵の名無しさん
20/03/04 20:28:09 u9i0Khrz.net
>>983
受験を控えている高校生です
他のスレでアナタはコンプ性Fラン症候群という不治の病であるとお聞きしました
高学歴の私と違ってFランであるアナタは今はどのようなお仕事をされてるんですかぁ?
1047:卵の名無しさん
20/03/04 22:38:01 gP3NYsfY.net
>>984
>973に答えて高学歴の学力をみせつけてみ!
1048:卵の名無しさん
20/04/24 23:02:54 CaIpZWAk.net
医療大麻 CBDオイル 匿名ユーザーレビュー特集
・なんで効くのかわからないけど大麻すごい
父の病気が、治りました。
医者が驚いてます。僕も驚いてます。父は喜んでます。母は得意げです。
↓↓
謎に迫る! URLリンク(plaza.rakuten.co.jp)
1049:卵の名無しさん
20/08/28 15:11:50 I97SIklY.net
num lockキーをback space(BS)に設定する。
Windows Registry Editor Version 5.00
[HKEY_LOCAL_MACHINE\SYSTEM\CurrentControlSet\Control\Keyboard Layout]
"Scancode Map"=hex:00,00,00,00,00,00,00,00,02,00,00,00,0e,00,45,00,00,00,00,00
URLリンク(www.souichi.club)
1050:卵の名無しさん
20/09/28 07:18:00.74 nQ3vcs5K.net
- ショートカットを表示しない
[HKEY_CURRENT_USER\Software\Microsoft\Windows\CurrentVersion\Explorer]
"link"=hex:00,00,00,00
1051:卵の名無しさん
20/09/28 07:18:23.22 nQ3vcs5K.net
- ショートカットを表示しない
[HKEY_CURRENT_USER\Software\Microsoft\Windows\CurrentVersion\Explorer]
"link"=hex:00,00,00,00
1052:卵の名無しさん
20/12/27 22:55:28.44 Wj30LPWM.net
水の比重1、重力加速度9.8m/(s^2)とすると
1N=1/9.8kg重
1hPa=100N/m^2
100/9.8kg重/(100cm)^2
100/9.8=10.19716
100/9.8(kg重)/(100)^2(cm^2)
100/9.8*1000(g重)/100^2(cm^2)
100/9.8*1000/100^2(g重/cm^2)
1053:卵の名無しさん
20/12/29 00:37:18.34 RkgAFAj8.net
うりゅうひろゆきの本スレにようこそ!
このウンコスレは基本的にさらしアゲて使用します!
ココは医師真似事務員ジイさんが、医師妬みの挙句、テメェで勝手に自称で医科歯科大を卒業した医師って設定で
妄想願望日記を発表するスレです
お医者さんのマトモなレスなど全く不要!
ゴミ箱まがいの公衆便所掃き溜め痰壷下品スレとして使用しましょうね
そんなわけで、いまだ懲りない恥さらし(自称)医科歯科事務員ジジイ登場~~
│ 彡川川川三三三ミ~ プウゥ~ン
| 川|川/ \|~ ポワ~ン ________
| ∥|∥ ◎---◎|~ /
| 川川∥ 3 ヽ~ < 僕はうりゅうひろゆき
| 川川 ∴)д(∴)~ \________
| 川川 ~ /~ カタカタカタ
| 川川∥ ~ /∥ _____
| 川川川
1054:川___/∥ | | ̄ ̄\ \ | / \__| | | ̄ ̄| | / \医科歯科 | | |__| | | \ |つ |__|__/ / | / 医師うらやましいぜ・・ | ̄ ̄ ̄ ̄| 〔 ̄ ̄〕 よし、医科歯科の医者と称して 私立のお医者様を攻撃だ!なんたって、医科歯科の看板借りてんだからな え?証明???自称だから医科歯科の学生証、卒業証書あるわけねえよ 同期の医師あげてみろ、ったて、そもそも医師じゃねえし居ねえよバカ でも病院でもらった、医科歯科の封筒はもってま~す!これを証拠にしてやるぜ 専門医?欲しいけど・・いいや!そんなの必要ねえ! 俺は事務員だし
1055:卵の名無しさん
21/01/12 13:57:21.80 5eSow97R.net
ウリュ爺見てるか?
1056:卵の名無しさん
21/01/17 11:43:16.68 0+l4A2WC.net
r1=100
n1=10000
r2=500
n2=10000
k=1e6
va=rbeta(k,1+r1,1+n1-r1) # vacccine
pl=rbeta(k,1+r2,1+n2-r2) # placebo
rr=pl/va
summary(rr)
quantile(rr,c(0.025,0.975))
NNT=1/(pl - va)
summary(NNT)
quantile(NNT,c(0.025,0.975))
1057:卵の名無しさん
21/01/17 14:22:04.75 0+l4A2WC.net
20000人の被験者の半数にはワクチンを、半数には偽薬を
注射し、1ヶ月後に調べたところ、ワクチン群には10人、
偽薬群には50人の感染者が出た。
このワクチン接種により感染確率が1/4未満になる確率はいくらか?
r1=10
n1=1000
r2=50
n2=1000
k=1e6
va=rbeta(k,1+r1,1+n1-r1) # vacccine
pl=rbeta(k,1+r2,1+n2-r2) # placebo
rr=pl/va
summary(rr)
quantile(rr,c(0.025,0.975))
BEST::plotPost(rr,xlab='risk ratio')
HDInterval::hdi(rr)
mean(1/rr < 1/4)
NNT=1/(pl - va)
summary(NNT)
quantile(NNT,c(0.025,0.975))
BEST::plotPost(NNT,xlab='NNT',col=2)
HDInterval::hdi(NNT)
1058:卵の名無しさん
21/01/17 22:21:32.82 Phpvqw96.net
'
感染有 感染無
既往有 44 6570
既往無 409 13764
'
covid=
cbind(
replicate(44, c(1,1)),
replicate(409, c(1,0)),
replicate(6570, c(0,1)),
replicate(13764, c(0,0)))
COVID=t(covid)
colnames(COVID)=c('infection','history')
COVID19=as.data.frame(COVID)
glm=glm(infection~.,data=COVID19, family=binomial)
adj.or=exp(glm$coef)
round(cbind(adj.or, exp(confint(glm))),3)
k=1e6
neg=rbeta(k,1+r1,1+n1-r1) # negative history
pos=rbeta(k,1+r2,1+n2-r2) # positive history
neg.o=neg/(1-neg)
pos.o=pos/(1-pos)
or=pos.o/neg.o
quantile(or,c(0.025,0.975))
summary(or)
1059:卵の名無しさん
21/01/18 23:51:49.23 cEjl/X0L.net
医者になる前はsoftware engineerだった俺には
RなんぞVBレベル痴呆言語にしか思えない。
そもそもdata.frameとか取ってつけたような拡張が
本当にうぜえ。S言語自体が80年代の設計だから
マトリックス操作がcbindやらなんちゃらbindやら
オプション満載の超原始的。
んな言語恥ずかしいレベルやな。
1060:卵の名無しさん
21/01/19 21:49:17.29 kd8IADuM.net
ウリュウ=プログラムおじさん
社会にも5chでも除け者にされ居場所はない。
1061:卵の名無しさん
21/01/19 21:49:50.10 kd8IADuM.net
ウリュウに待つのは医者コンプに塗れながらの孤独死である。
1062:卵の名無しさん
21/01/19 21:50:52.32 kd8IADuM.net
東医の佐野君→若くて実家がお金持ちで将来は医者
ウリュウ→金なし、医師免許なし、将来なし
1063:卵の名無しさん
21/01/19 21:51:27.13 kd8IADuM.net
私立医コンプレックスのウリュウのジジイにつける薬など
1064:ない。
1065:1001
Over 1000 Thread.net
このスレッドは1000を超えました。
新しいスレッドを立ててください。
life time: 811日 23時間 32分 1秒
1066:過去ログ ★
[過去ログ]
■ このスレッドは過去ログ倉庫に格納されています