風俗嬢ですが再起して頑張り医学部に合格しました [転載禁止]©2ch.net
■ このスレッドは過去ログ倉庫に格納されています
風俗板にて
風俗辞めた後に医学部に合格したと嘘をついて袋叩きにされてる様子
http://aoi.bbspink.com/test/read.cgi/club/1424999821/724
760 ：名無しさん＠ピンキー：2015/04/03(金) 17:34:23.58 ID:hEZXryQ10
>>724
アホか、お前
＞医者になった嬢も居る、
そんなの絶対いないわ、ボケ。
お前は医者と看護師の区別もつかんのか(大笑) この誤った英作文Your thing〜は入試の学力考査が適正に行われていない証であり、
誤っていることすらわからないのは、国試に関係ない分野での大学教育も健全に行われていないことの反映である。
http://imagizer.imageshack.com/img924/2411/sYBNlI.jpg
医大入学を金で買う をベイズ推計してみよう。
帰無仮説H0:学力考査は公正に行われている
対立仮設H1:縁故加点など不正が行われている
H0の成り立つ事前確率をP0
H1の成り立つ事前確率をP1
p : 学力考査が公正であった場合にYour thing英作文をする確率P(Y|H0)
q : 学力考査が不正であった場合にYour thing英作文をする確率P(Y|H1)
P(H1|Y)=P(Y|H1)P(H1)/P(Y)=P(H1)P(Y|H1)/P(H0)(P(Y|H0)+P(H1)P(Y|H1) (ベイズの公式)から
H1の事後確率（裏口入試存在確率）は P1q/(P0p＋P1q) である。
情報がないという前提でP0=P1=0.5と仮定しよう。（縁故加点合格の証言もあるが）
公正な学力考査で選抜された学生がYour thing英作文をする確率を0.01としてみる（実際はもっと低いはず）
ド底辺特殊シリツ医大卒にYour thing英作文をする底抜けバカの占める割合を存じ上げないのでグラフにしてみた。
http://i.imgur.com/5vesmwb.jpg # OR = ad/bc
# q1 = a/(a+c) 疾病群での曝露リスク
# q0 = b/(b+d) 対照群での曝露リスク
# RR = q1/q0
# RR =a/(a+c) ÷ b/(b+d)=a(b+d)÷b(a+c)=(ab+ad)/(ab+bc)
# =(ab+bc*OR)/(ab+bc) ∵ ad=bc*OR
# =(c*OR+a)/(a+c)
# =OR*c/(a+c)+a/(a+c)
# =OR*( 1-a/(a+c) ) + a/(a+c)
# =OR*(1-q1) + q1
# RR = a/(a+c) ÷ b/(b+d)=a(b+d)÷b(a+c)=(ab+ad)/(ab+bc)
# = (ab+ad)/(ab+ad/OR) ∵ bc=ad/OR
# = (b+d) /(b + d/OR) # (b+d) で割る
# =1 / {b/(b+d) + d/[(b+d)/OR]}
# =OR/{OR*b/(b+d) + d/(b+d)} # ORを掛ける
# =OR/{OR*b/(b+d) + (b+d)/(b+d)-b/(b+d)}
# =OR/(OR*q0 + 1-q0) or2p=function(a,b,c,d,x){
OR=(a/b)/(c/d)
SE=sqrt(1/a+1/b+1/c+1/d)
pu=(1-pnorm(log(x/OR)/SE))*2
pl=(1-pnorm(-log(x/OR)/SE))*2
p.value=ifelse(pl<=1,pl,pu)
return(p.value)
}
par(mfrow=c(2,1))
# 男子学生
curve(or2p(206,179,430,86,x),from=0.01,to=2,lwd=2,log="x",xlab="OR",y="p.value",las=1)
or2p(206,179,430,86,1)
abline(h=0.05,lty=3,col=2)
abline(v=1,lty=2,col=4)
legend("topleft",bty="n",legend=c("p=0.05","OR=1"),lty=c(3,2),col=c(2,4))
Epi::twoby2(dat.m)
chisq.test(dat.m,correct=FALSE)
fisher.test(dat.m)
# 女子学生
curve(or2p(26,13,317,47,x),from=0.01,to=2,lwd=2,log="x",xlab="OR",y="p.value",las=1)
or2p(26,13,317,47,1)
abline(h=0.05,lty=3,col=2)
abline(v=1,lty=2,col=4)
legend("topleft",bty="n",legend=c("p=0.05","OR=1"),lty=c(3,2),col=c(2,4))
dat.f=matrix(c(26,13,317,47),2)
Epi::twoby2(dat.f)
chisq.test(dat.f,correct=FALSE)
fisher.test(dat.f) xchg=function(a,b,bit=1){ # a, bをランダムにbitずつランダム交換
flg=runif(1)>0.5
if(flg){
chg=ifelse(b>=bit,bit,b)
return(chg)
}else{
chg=ifelse(a>=bit,-bit,-a)
return(chg)
}
}
rchg=function(x,bit=1){
n=length(x)
idx=sample(n,2)
chg=xchg(x[idx[1]],x[idx[2]],bit)
x[idx[1]]=x[idx[1]]+chg
x[idx[2]]=x[idx[2]]-chg
return(x)
}
mono=c(rep(0,99),500)
hist(mono,freq=FALSE)
for(i in 1:10000) mono=rchg(exp)
hist(mono,freq=FALSE)
library(fitdistrplus)
fit=fitdist(mono,"exp") ; fit$est
curve(dexp(x,rate=fit$est),add=TRUE) HDI<-function(x,conf.level=0.95){
n=length(x)
xs=x/sum(x)
lim=1-conf.level
cul=0 ; i=1 ; j=n
while(cul<lim){
if(xs[i]<xs[j]){
cul=cul+xs[i]
i=i+1
} else {
cul=cul+xs[j]
j=j-1
}
}
ans=c(i,j)/n
return(ans)
} ##
dx <- density(x)
dn <- cumsum(dx$y)/sum(dx$y)
li <- which(dn>=0.025)[1]
ui <- which(dn>=0.975)[1]
xx[li] ; xx[ui]
dx$x[c(li,ui)]
##
HDI<-function(x,conf.level=0.95){
n=length(x)
xs=x/sum(x)
lim=1-conf.level
cul=0 ; i=1 ; j=n
while(cul<lim){
if(xs[i]<xs[j]){
cul=cul+xs[i]
i=i+1
} else {
cul=cul+xs[j]
j=j-1
}
}
ans=c(i,j)/n
return(ans)
} 汚らわしい 作品ばかりで 恥かしい 裏口バカが 題材ゆえに
ド底辺 出身校すら 名のれない 裏口入学 恥ずかしいから？
まやかしの 学力考査の ド底辺 書いた英語は Your thing.
いつの世も 特殊民の ド底辺 受験するのは バカ子息のみ
補欠です 寄付金いくら 払えます？ 事前予告は お上も容認。
本人も 認めているよ 私立馬鹿 裏口バカと 同義なり
http://i.imgur.com/rDRhIsg.jpg
法外な カネを払って ド底辺 裏口バカと 蔑まれるに
世の中で 最も恥ずべき 資格かな 屈辱ゆえに 校名名のれず
世の中の 誰もが知ってる ド底辺 学力不足を 金で補う
Your thing 一目瞭然 ド底辺 高校入試も 受からんぞ nbinom.test=function(r,n,p=0.5,conf.level=0.95){
HDI<-function(x,cl=conf.level){
n=length(x)
xs=x/sum(x)
lim=1-cl
cul=0 ; i=1 ; j=n
while(cul<lim){
if(xs[i]<xs[j]){
cul=cul+xs[i]
i=i+1
} else {
cul=cul+xs[j]
j=j-1
}
}
ans=c(i,j)/n
return(ans)
}
L=function(x) dnbinom(n-r,r,x)
sL=function(x) L(x)/integrate(L,0,1)$value
pL=function(x) integrate(sL,0,x)$value
xx=seq(0,1,le=10001) ; xx=xx[-1]
yy=sapply(xx,sL)
hdi=HDI(yy,cl=conf.level)
fu0=function(x,u0) pL(x)-u0
ci=c(lower=uniroot(fu0,u0=(1-conf.level)/2,c(1e-16,1))$root,
upper=uniroot(fu0,u0=(1-(1-conf.level)/2),c(1e-16,1))$root)
p.value=1-pnbinom(n-1-r,size=r,p)
round(data.frame(ci,hdi,p.value),4)
} nbinom.test=function(r,n,p=0.5,conf.level=0.95){
HDI<-function(x,cl=conf.level){
n=length(x)
xs=x/sum(x)
lim=1-cl
cul=0 ; i=1 ; j=n
while(cul<lim){
if(xs[i]<xs[j]){
cul=cul+xs[i]
i=i+1
} else {
cul=cul+xs[j]
j=j-1
}
}
ans=c(i,j)/n
return(ans)
}
L=function(x) dnbinom(n-r,r,x)
sL=function(x) L(x)/integrate(L,0,1)$value
pL=function(x) integrate(sL,0,x)$value
xx=seq(0,1,le=10001) ; xx=xx[-1]
yy=sapply(xx,sL)
plot(xx,yy,type="l",xlab="��",ylab="Density",lwd=2,las=1, frame.plot=FALSE)
hdi=HDI(yy,cl=conf.level)
fu0=function(x,u0) pL(x)-u0
ci=c(lower=uniroot(fu0,u0=(1-conf.level)/2,c(1e-16,1))$root,
upper=uniroot(fu0,u0=(1-(1-conf.level)/2),c(1e-16,1))$root)
p.value=1-pnbinom(n-1-r,size=r,p)
round(data.frame(ci,hdi,p.value),4)
} nbinom.ci=function(r,n,p=0.5,conf.level=0.95,...){
HDI<-function(x,cl=conf.level){
n=length(x)
xs=x/sum(x)
lim=1-cl
cul=0 ; i=1 ; j=n
while(cul<lim){
if(xs[i]<xs[j]){
cul=cul+xs[i]
i=i+1
} else {
cul=cul+xs[j]
j=j-1
}
}
ans=c(i,j)/n
return(ans)
}
L=function(x) dnbinom(n-r,r,x)
sL=function(x) L(x)/integrate(L,0,1)$value
xx=seq(0,1,le=1001) ; xx=xx[-1]
yy=sapply(xx,sL)
plot(xx,yy,type="l",xlab="θ",ylab="Density",lwd=2,las=1, frame.plot=FALSE,...)
hdi=HDI(yy,cl=conf.level)
psL=cumsum(yy)/sum(yy)
ci=c(lower=which(psL>(1-conf.level)/2)[1]/length(yy),
upper=which(psL>(1-(1-conf.level)/2))[1]/length(yy))
p.value=1-pnbinom(n-1-r,size=r,p)
round(data.frame(ci,hdi,p.value),4)
} bpt=function(t) sqrt(m/(2*pi*a^2*t^3))*exp(-(t-m)^2/(2*a^2*m*t))
x=c(11*12+5,12*12+5,15*12+5,15*12+9)
x=sort(x) ; x
dx=diff(x) ; dx
m=mean(dx) ; m
a=sd(dx)/mean(dx) ; a
bpt=function(t) sqrt(m/(2*pi*a^2*t^3))*exp(-(t-m)^2/(2*a^2*m*t))
tt=seq(0,10*12,by=1)
y=sapply(tt,bpt)
dev.off()
plot(tt,y,type="l",lwd=1, xlab="months since May 2011",ylab="density")
17*12-(11*12+5) # 67
16*12-(11*12+5) # 55
integrate(bpt,55,67)
xxx=seq(55,120,le=480)
segments(xxx,0,xxx,sapply(xxx,bpt),col="tomato")
xx=seq(55,67,le=24)
segments(xx,0,xx,sapply(xx,bpt))
integrate(bpt,55,67)$value/integrate(bpt,55,Inf)$value p-hackingのネタ
女子大生100人と女子高生100人の胸囲を測定して
前者が平均82 , 標準偏差3
後者が平均81 , 標準偏差3
であったとする。
統計的に有意差があるといえるか？
## t.test & Welch's t.test
ttest = function(n1,m1,s1, n2,m2,s2){ # ni:標本数 mi:平均 si:標準偏差
n=n1+n2-2
u=((n1-1)*s1^2+(n2-1)*s2^2)/n
t=(m1-m2)/sqrt(u/n1+u/n2)
r=cat("Equal Variance:\n\t",sep="")
r=cat(r,"t = ",round(t,2),",\tdf = ", n ,",\tp = ", 2*pt(-abs(t),n), "\n" , sep="")
t=(m1-m2)/sqrt(s1^2/n1 + s2^2/n2)
n=(s1^2/n1+s2^2/n2)^2 / ( (s1^2/n1)^2/(n1-1) + (s2^2/n2)^2/(n2-1))
r=cat(r,"Unequal variance:\n\t",sep="")
cat(r,"t = ",round(t,2),",\tdf = ", n ,",\tp = ", 2*pt(-abs(t),n), "\n" , sep="")
}
n1=100 ; m1=82 ; s1=3 n2=100 ; m2=81 ; s2=3
x1=scale(rnorm(n1))*s1+m1
x2=scale(rnorm(n2))*s2+m2
t.test(x1,x2)
t.test(x1,x2,var=TRUE)
ttest(n1,m1,s1, n2,m2,s2) >>655
, -――――――‐- ､
, ｀丶､
／ ＼
／ ＼
/ ｌ ヽ
ﾊ l ､ ｌ ｌ .
| .ｌ / 、 | ＼ |、 | | ＿|＿ ＿|＿ ＿＿＿ -―フ
| /| ./ ヽ | ＼ ｜＼ | | ノ ＼ | ＼ | ＼ ／ ∠
| / |/ ./ ＼ | ｀ ーヽ‐ ､ ヽ | Ｘ | ｜ | | / ）
ﾉ /, ― ' ＼| , =＝＝ミ. V ､ ／ ＼ | ｜ ＼／ ＼ ○ノ
/ { / , =＝ミ､ 〃{//////ﾊヽ | ､
/ V 〃////ﾊ V/////:ﾉ | ｌへ ヽ
/／ ﾍ { {/////ﾉ =＝＝= | |＾ } ＼ ＿|＿ ＿|＿ ＿＿＿ -―フ
| ､ ヾ =＝ / / | ' ./ ﾄー' ノ ＼ | ＼ | ＼ ／ ∠
| } / / ハ | ./／ | 、 | Ｘ | ｜ | | / ）
. ﾊ { ‐r―‐ ´ } ／|/ .| ＼| ／ ＼ | ｜ ＼／ ＼ ○ノ
∨ ヽ、 ､＿＿ ノ , イ ﾊ. |
| ﾊ ﾊ＞‐/^^^^ヽ‐r−r:＜ﾊ:ミ＜Λ/∨ ヽ|
|/ ∨ ∨/＝＝=/ヽ_><_／ | ＼
/ ./> | ∨ ./ ＜ / | 理系男「”eのx乗を不定積分”しようぜ！」
理系女「何のこと？」
理系男「∫eXdx…。これのdxより前のところ見てごらん。何かに見えないかい？」
理系女「sexよね？」
理系男「気持ち良くなりたいだろ？」
理系女「ええ。」
理系男「じゃあ実際に積分してみるんだ。さあ、はやく」
理系女「eX+Cになったわ」
理系男「どうだ、気持ちよくなっただろ？eX+Cの「+」はプラスじゃなくて「足す」と読むんだ」
理系女「…エクスタシー！！」
理系男「どうだ、気持ちよくなっただろ？」
理系男女「イク-----ッ！！」 count=function(N=(1:10)*10,.pch=1:20,.col=1:8,...){
n=sample(N,1)
x=rnorm(n)
y=rnorm(n)
plot(x,y,axes=FALSE,ann=FALSE,pch=.pch,col=.col,...)
n
}
par(mfrow=c(2,2))
count(,19,1)
count(,19,1,cex=1:2)
count(,19,)
count(,19,cex=1:2) 当直医のスレ Part.17 [無断転載禁止]©2ch.net
http://potato.2ch.net/test/read.cgi/hosp/1480133830/817
817 名前：卵の名無しさん[] 投稿日：2016/12/17(土) 19:17:54.68 ID:3Oft5frV
>>786
彼の目的が分らない。専門医もとらず学位もとらず、
管理職にも就かせてもらえず、医局を首になったか何かのトラウマ？
それで一生懸命に私学を叩いては喜んでる。
異常としか言いようがない。
俺はあんたが異常だと思うわ。
管理職なんか成りたいか？
所詮、病院の経営者から雇われてるだけの雇われもんじゃんw
雇われ院長になって安い給料をもらってまわりからちやほやされて
喜んでる医者の方が異常だろ。器が小さすぎだろ（笑）。
俺は内視鏡バイトや麻酔のバイトやって金を稼ぎまくってる自称医科歯科大学卒の
爺さんに好感を持つわ。この爺さんの文章を読んでると結構頭良いし、もし内視鏡
や麻酔科のバイトができる技術を持ってると言うのが本当なら、相当に実力のある
医者だと思うぞ。 >>817
卒後就職したのは、麻酔、内視鏡、アンギオは外科医なら誰でもやるのが当然の地方の病院だった。
製品化される前のエトキシスクレロールで硬化療法とか消化器全般が外科の仕事だった。腋窩切開で気胸の手術くらいは自力でやっていた。
整形外科不在時は脱臼整復やギプス固定までせざるをえなかった。
帝王切開の助手は外科医の仕事だった。仮死状態で生まれた新生児の挿管までした
(暴れる元気はないけど首がすわってないから固定をしっかりすれば難しくない）。
今はリスクの少ない手技しかしなくなったというより、しなくてよい職場を選んでいる。
ESDで大腸穿孔の経験もある。後腹膜側だったから保存的に治癒。最近はEMRまでしかしない。 当直医のスレ Part.17 [無断転載禁止]©2ch.net
http://potato.2ch.net/test/read.cgi/hosp/1480133830/827
827 名前：卵の名無しさん[] 投稿日：2016/12/18(日) 01:59:11.53 ID:Yv/ypoQO
>>821
60歳代のウチの外科専門医を持ってる院長とは大違いだなw
ESDは言うに及ばずEMRどころかCSすら深部到達できないw
いつも大口叩いてるくせに救急処置が苦手で挿管も超下手くそで
CVすら取れないと知ってからは気の毒なのであまり責めないようにしているw
いつも専門医自慢してウザかったけど挿管が超下手くそでCV取れないのが
看護師らにバレてからは大人しくなったw pay=function(income101_172,u0=0)0.5916*income101_172+14884-u0
f=function(x)pay(x*10000/12)*12
f(1772) f4=function(N=40,n1=20,n2=20,.replace=TRUE){
x1=sample(N,n1,replace=.replace)
x2=sample(N,n2,replace=.replace)
xx12=rank(c(x1,x2))
xx1=xx12[1:n1]
xx2=xx12[(n1+1):(n1+n2)]
c(exactRankTests::wilcox.exact(x1,x2)$p.value,lawstat::brunner.munzel.test(x1,x2)$p.value,
t.test(x1,x2,var.equal=FALSE)$p.value,t.test(xx1,xx2)$p.value,wilcox.test(x1,x2)$p.value)
}
f4()
f4(.replace=FALSE)
r=replicate(1000,f4(N=100,n1=20,n2=20,.replace=FALSE))
dev.off()
par(mfrow=c(2,2))
plot(r[1,],r[2,],asp=1,xlab="Wilcox.exact",ylab="Brunner Munzel")
plot(r[4,],r[2,],asp=1,xlab="順位t.test(Welch)",ylab="Brunner Munzel")
plot(r[4,],r[3,],asp=1,xlab="順位t.test(Welch)",ylab="t.test(Welch)")
plot(r[2,],r[3,],asp=1,xlab="Brunner Munzel",ylab="t.test(Welch)")
plot(r[1,],r[5,],asp=1,xlab="wilcox.exact",ylab="wilcox")
r=replicate(1000,f4(N=5,n1=20,n2=20,.replace=TRUE))
dev.off()
par(mfrow=c(2,2))
plot(r[1,],r[2,],asp=1,xlab="Wilcox.exact",ylab="Brunner Munzel")
plot(r[4,],r[2,],asp=1,xlab="順位t.test(Welch)",ylab="Brunner Munzel")
plot(r[4,],r[3,],asp=1,xlab="順位t.test(Welch)",ylab="t.test(Welch)")
plot(r[2,],r[3,],asp=1,xlab="Brunner Munzel",ylab="t.test(Welch)")
plot(r[1,],r[5,],asp=1,xlab="wilcox.exact",ylab="wilcox") 医師もMRも幸せにする患者のための情報吟味: ディオバン事件以降の臨床研究リテラシー 臨床医なんて接客業。
客がいなけりゃただの人 愛想良く媚びることを旨とした客商売なので
役者・芸者と同じ。
ド底辺特殊シリツ医大卒でも客商売がうまければ稼げる。
ド底辺特殊シリツ医大卒の悲哀はこれ、
＞＞
私は昭和の時代に大学受験したけど、昔は今よりも差別感が凄く、慶応以外の私立医は特殊民のための特殊学校というイメージで開業医のバカ息子以外は誰も受験しようとすらしなかった。
常識的に考えて、数千万という法外な金を払って、しかも同業者からも患者からもバカだの裏口だのと散々罵られるのをわかって好き好んで私立医に行く同級生は一人もいませんでした。
本人には面と向かっては言わないけれど、俺くらいの年代の人間は、おそらくは８−９割は私立卒を今でも「何偉そうなこと抜かしてるんだ、この裏口バカが」と心の底で軽蔑し、嘲笑しているよ。当の本人には面と向かっては絶対にそんなことは言わないけどね。
＜＜ ## y1,y2が独立にポアソン分布に従うとき
# y1 ~ pois(λ1)
# y2 ~ pois(λ2)
# N = y1 + y2
# y1 ~ Binomial(N,λ1/(λ1+λ2))
n=100
lambda=sample(0:n,2)
l1=lambda[1] ; l2=lambda[2]
y1=rpois(1000,l1)
y2=rpois(1000,l2)
N=y1+y2
yy1=rbinom(1000,N,l1/(l1+l2))
plot(density(y1),lwd=2)
lines(density(yy1),lwd=2,col=4)
hist(y1,freq=FALSE)
lines(density(yy1),lwd=2,col=4)
http://i.imgur.com/y2kXM6r.jpg ## bernoulli.stan
data{
int N;
int<lower=0,upper=1> Y[N];
}
parameters{
real<lower=0,upper=1> p;
}
model{
Y ~ bernoulli(p);
}
##
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
stanmodel <- stan_model('bernoulli.stan')
Sn <- function(S,n,...){
dat = sample(c(rep(1,S),rep(0,n-S)))
d = data.frame(Y=dat)
data=list(N=nrow(d),Y=d$Y)
fit=rstan::sampling(stanmodel,data=data)
print(fit,pars="p",...)
}
Sn(3,10)
Sn(30,100)
Sn(5,100) library(PSAgraphics)
data('lindner')
library(Matching)
head(lindner)
summary(glm(lifepres ~ .-1, data=lindner))
Y=lindner$lifepres
Tre=lindner$abcix
ps=glm(abcix ~ .,data=lindner[,-1], family=binomial(link="logit"))$fitted
mout.lindner=Match(Y,Tre,ps)
summary(mout.lindner)
ivec1 <- lindner$abcix
ivec2 <- rep(1,nrow(lindner)) - ivec1
ivec <- cbind(ivec1, ivec2)
iestp1 <- (ivec1/ps)/mean(ivec1)
iestp2 <- (ivec2/ps)/mean(ivec2)
iestp <- iestp1 + iestp2
ipwe <- lm(Y ~ ivec -1, weights=iestp, data=lindner)
summary(ipwe)
library(survey)
psw <- svydesign(ids = ~ 1 , weights=iestp, data=lindner)
ipwsu <- svyglm(lifepres ~ abcix, design = psw)
summary(ipwsu) p=0.4 # prevalence
se=0.62 # sensitivity
sp=0.98 # specificity
a=p/(1-p) # odds
b=a*se/(1-sp) # *LH(+)
b/(1+b) # ppv
c=a*(1-se)/sp # *LH(-)
1 - c/(1+c) # npv ## p.351 TABLE19-3
A1=45
A0=94
B1=257
B0=945
dat=matrix(c(A1,B1,A0,B0),2)
colnames(dat)=c('X1','X0')
rownames(dat)=c('D1','D0')
dat
Pz1=c(0.40,0.55,0.70,0.45,0.60,0.75) #暴露喫煙率
Pz0=c(0.30,0.45,0.60,0.25,0.40,0.55) #非暴露喫煙率
ORxz = (Pz1/(1-Pz1))/(Pz0/(1-Pz0)) ; round(ORxz,2)
ORdz=c(5,10,15) # 病気の有無による喫煙オッズ比(暴露によらず一致)
n=length(ORdz)
m=length(ORxz)
B11=Pz1*B1 ; B11 #無病暴露喫煙者数
B01=Pz0*B0 ; B01 #無病非暴露喫煙者数
A11=A01=matrix(numeric(n*m),ncol=n)
ORdx=numeric(n)
for(i in 1:n){
A11[,i] = ORdz[i]*A1*B11/(ORdz[i]*B11 + B1-B11) #有病暴露喫煙者数
A01[,i] = ORdz[i]*A0*B01/(ORdz[i]*B01 + B0-B01) #有病非暴露喫煙者数
ORdx = (A11*B01)/(A01*B11)
# ORdxz0 = (A1-A11)*(B0-B01)/((A0-A01)*(B1-B11))
}
round(ORdx,2) # ある文字列が迷惑メールに含まれる確率をp1[i],通常メールに含まれる確率をp0[i]とする。
x1=0.1 # 初期値
p1=rep(0.1,7)
p0=rep(0.01,7)
L.J.Savage <- function(x1){
Mail <- function(x,i){
p1[i]*x/(p1[i]*x +p0[i]*(1-x))
}
n=length(p1)
X=numeric(n)
X[1]=x1
for(i in 1:n){
X[i+1]=Mail(X[i],i)
}
return(X)
}
L.J.Savage(0.1)
xx=c(0.01,0.05,0.20,0.30,0.40,0.50)
sapply(xx,L.J.Savage) >>692
当直スレのレボジトリに使ってんの。
直接スクリプトを貼ると読めない人には迷惑だからね。 ##
# あるウイルス疾患の診断に迅速キット、発熱、関節痛の
# 感度は80% 90%,70%、特異度は99%,50%,60%とする。
# 迅速キット：陰性、発熱：陽性、関節痛：陽性であった。
# 診察前確率を50%としたとき、個々の検査陽性確率は独立という前提で診察後確率を求めよ。
x=0.5 # prevalene
sn=c(0.8,0.9,0.7) # sensitivity
sp=c(0.99,0.5,0.6) # specificity
E=c(0,1,1) # exam result
TP=sn
FP=1-sp
pLH=TP/FP
nLH=(1-TP)/(1-FP)
p2o <- function(p) p/(1-p)
o2p <- function(o) o/(1+o)
LH=prod(ifelse(E,pLH,nLH))
o2p(p2o(x)*LH) このあたりの話はMcGeeに書いてあるはず。
俺は理解できているから尤度比の参照にしか使ってない。
McGeeのデータは尤度比の記載のみで
感度特異度は割愛されているが次の式で算出できる。
p:陽性尤度比 q:陰性尤度比
感度=p(1-q)/(p-q)
特異度=(p-1)/(p-q) # 3項目のうち少なくとも２項目が陽性で疾患ありとしたときの感度
n=length(sn)
sn2=numeric(n)
for(i in 1:n){
sn2[i] <- (1-sn[i])*prod(sn[-i])
}
sum(sn2)+prod(sn)
# 3項目のうち少なくとも２項目が陽性で疾患ありとしたときの特異度
n=length(sp)
sp2=numeric(n)
sp2
for(i in 1:n){
sp2[i] <- (1-sp[i])*prod(sp[-i])
}
sum(sp2)+prod(sp)
# 3項目がすべて陽性で疾患ありとしたときの特異度
1-prod(1-sp) ##
# Aa,Ba,Ca : 恩赦事象 1/3
# Bd,Cd : 死刑宣告 1/2
# p(Aa|Bd) = p(Bd|Aa)*p(Aa) / p(Bd)
# where p(Bd) = p(Bd|Aa) *P(Aa) + p(Bd|Ba)*p(Ba) + p(Bd|Ca)*p(Ca) = 0.5*0.33 + 0*0.33 + 1*0.33 =1/2*1/3 + 0 + 1*1/3 = 3/6 = 1/2
# ∴ p(Aa|Bd) = (1/6)/(1/2) = 1/3 ##
# Aa,Ba,Ca : 恩赦事象 1/3
# Bd,Cd : 死刑宣告 1/2
# p(Aa|Bd) = p(Bd|Aa)*p(Aa) / p(Bd)
# where p(Bd) = p(Bd|Aa)*P(Aa) + p(Bd|Ba)*p(Ba) + p(Bd|Ca)*p(Ca)
= (1/2) *(1/3) + 0*(1/3) + 1*1/3
= 1/2
# ∴ p(Aa|Bd) = p(Bd|Aa)*p(Aa) / p(Bd)
= ((1/2)*(1/3)) / (1/2) = 1/3 当然ながら、迷惑メールフィルターは非迷惑メールの選択にも使える。
このスレのある投稿に含まれていた3つのキーワードは次の確率で国立卒の投稿（国立投稿）、ド底辺特殊シリツ医大卒の投稿に含まれるとする。
キーワード 国立卒 ド底辺特殊シリツ医大卒
裏口 0.1 0.001
底辺 0.3 0.01
ジジイ 0.001 0.25
この3つのキーワードを含む投稿がド底辺特殊シリツ医大卒の投稿である確率を述べよ。
事前確率は理由不在の原則に従って0.5とする。 # p(Aa|Bd) = p(Bd|Aa)*p(Aa) / p(Bd)
# where p(Bd) = p(Bd|Aa)*P(Aa) + p(Bd|Ba)*p(Ba) + p(Bd|Ca)*p(Ca)
# = (1/2) *(1/3) + 0*(1/3) + 1*1/3
# = 1/2
# ∴ p(Aa|Bd) = p(Bd|Aa)*p(Aa) / p(Bd)
# = ((1/2)*(1/3)) / (1/2) = 1/3
# P(Ca|Bd) = P(Bd|Ca)*P(Ca)/P(Bd)
# = ( 1 * (1/3))/ (1/2) = 2/3
N=10^5
ABC=c('A','B','C')
BC=c('B','C')
Amnesty=numeric(N)
Declare=numeric(N)
for(i in 1:N){
Amnesty[i]=sample(ABC,1)
switch(Amnesty[i],
'A'=Declare[i] <- sample(BC,1),
'B'=Declare[i] <- 'C',
'C'=Declare[i] <- 'B'
)
}
mean(Amnesty[which(Declare=='B')]=='A') #p(Aa|Bd)
mean(Amnesty[which(Declare=='B')]=='B') #p(Ba|Bd)
mean(Amnesty[which(Declare=='B')]=='C') #p(Ca|Bd) >>701
問題に答える頭もないの？
こういう問題は臨床でも必要だろ？
あるウイルス疾患の診断に迅速キット、発熱、関節痛の
感度は80% 90%,70%、特異度は99%,50%,60%とする。
迅速キット：陰性、発熱：陽性、関節痛：陽性であった。
診察前確率を50%としたとき、個々の検査陽性確率は独立という前提で診察後確率を求めよ。 私立医大の偏差値
http://2chreport.net/rank_01.htm をみると
ド底辺特殊シリツ医大の公称の偏差値は最低でも57となっている。
成績が正規分布と仮定して偏差値57は上位24%に位置する。
そのような学力の学生がこのような英作文をして
http://imagizer.imageshack.com/img924/2411/sYBNlI.jpg
しかも、間違いを訂正できない確率は極めて０に近いはずである。
100人にひとりそのようなバカがいたとするとその偏差値は下位１％なので26.7となる。
偏差値57以上の学生がYour thing英作文をする確率を1000分の1(実際はもっと少ないであろう)と大きめに設定する。
ここでベイズの公式
P(裏口入学|誤英作文) = P(誤英作文|裏口入学)* P(裏口入学)/ { P（誤英作文|裏口入学）* P(裏口入学) + P(誤英作文|正規入学)*P(正規入学) }
P(正規入学) = 1 - P(裏口入学), P（誤英作文|裏口入学) = 1
理由不十分の原則(principal of insufficient reason) http://mathworld.wolfram.com/PrincipleofInsufficientReason.html
にしたがって
P(裏口入学)=P(正規入学)=0.5とする。
P(誤英作文|正規入学)=0.001の設定であるから
P(裏口入学|誤英作文) = 0.5/(0.5 +0.001*0.5) = 0.999
と算出できる。
ド底辺特殊シリツ医大の公称最低偏差値が57が正しいとすると
このスレのド底辺特殊シリツ医大卒が裏口入学である可能性は99.9％である。 Q.E.D. ■ このスレッドは過去ログ倉庫に格納されています
