臨床統計もおもしろいですよ、その1 [無断転載禁止]©2ch.net

1卵の名無しさん2017/05/03(水) 20:04:54.62ID:0YB5L7xG
 
 内科認定医受験の最低限の知識、
 製薬会社の示してくる臨床データ、
 論文の考察、
 論文を書くときの正当性、
 というのが、臨床統計の今までの目的の大きい部分でしたが、
 
 AI=機械学習の基本も、結局は統計学と確率に支配されます。
 そういう雑多な話をするスレです。
 

496卵の名無しさん2018/05/17(木) 07:09:30.35ID:1zxWbYbO
calc.FPR.p2 <- function(r1,r2,n1,n2,alpha=0.05){
p.val=prop.test(c(r1,r2),c(n1,n2))$p.value
k=length(c(r1,r2))
df=k-1
p1=r1/n1 ; p2=r2/n2
Pi=c(r1/n1,r2/n2)
theta.i=asin(sqrt(Pi))
delta=4*var(theta.i)*(k-1)
N=2/(1/n1+1/n2)
ncp1=N*delta
power.chi=pchisq(qchisq(1-alpha,df),df,ncp1,lower=FALSE)
qcrit=qchisq(max(p.val,1-p.val),df,ncp=0)
curve(dchisq(x,df),0,ncp1*3,xlab=quote(chi),ylab='Density',bty='n') # H0
curve(dchisq(x,df,ncp1),add=TRUE,lty=2,col=2) # H1
abline(v=qcrit,col='gray')
legend('top',bty='n',legend=c('H0','H1'),col=1:2,lty=1:2)
y0=dchisq(qcrit,df,0)
y1=dchisq(qcrit,df,ncp=ncp1)
FPR.equal=y0/(y0+y1)
FPR.less=p.val/(p.val+power.chi)
VAL=c(power=power.chi,p.value=p.val,FPR.equal=FPR.equal,FPR.less=FPR.less)
print(VAL,digits=3)
invisible(VAL)
}

calc.FPR.p2(95,85,100,110)

497卵の名無しさん2018/05/17(木) 10:16:58.64ID:Z8D8umCF
calc.FPR.p2 <- function(r1,r2,n1,n2,alpha=0.05){
p.val=prop.test(c(r1,r2),c(n1,n2))$p.value
k=length(c(r1,r2))
df=k-1
p1=r1/n1 ; p2=r2/n2
Pi=c(p1,p2)
theta.i=asin(sqrt(Pi)) # arcsine conversion
delta=4*var(theta.i)*(k-1) # sum of squares
N=2/(1/n1+1/n2) # harmonic mean, subcontrary mean
ncp1=N*delta # non-central parameter
power.chi=pchisq(qchisq(1-alpha,df),df,ncp1,lower=FALSE)
# qchisq(1-0.05,1): 3.841
qcrit=qchisq(max(p.val,1-p.val),df,ncp=0)
# qcrit = prop.test(c(r1,r2),c(n1,n2))$statistic
curve(dchisq(x,df),0,2*qcrit,xlab=quote(chi),ylab='Density',bty='n') # H0
curve(dchisq(x,df,ncp1),add=TRUE,lty=2,col=2) # H1
# power.chi=AUC of right half of the H1 curve
abline(v=qcrit,col='gray')
legend('top',bty='n',legend=c('H0','H1'),col=1:2,lty=1:2)
text(qcrit,0,round(qcrit,2))
y0=dchisq(qcrit,df,0)
y1=dchisq(qcrit,df,ncp=ncp1)
FPR.equal=y0/(y0+y1) # length ratio
FPR.less=p.val/(p.val+power.chi) # area ratio
FPR.alpha=alpha/(alpha+power.chi) # FPR before analysis
VAL=c(power=power.chi,p.value=p.val,FPR.equal=FPR.equal,
FPR.less=FPR.less,FPR.alpha=FPR.alpha)
print(VAL,digits=3)
invisible(VAL)
}
calc.FPR.p2(85,95,100,100, alpha=0.05)
calc.FPR.p2(85,95,100,100, alpha=0.01)

498卵の名無しさん2018/05/17(木) 11:17:41.13ID:Z8D8umCF
calc.FPR.p2 <- function(r1,r2,n1,n2,alpha=0.05){
p.val=prop.test(c(r1,r2),c(n1,n2))$p.value
k=length(c(r1,r2))
df=k-1
p1=r1/n1 ; p2=r2/n2
Pi=c(p1,p2)
theta.i=asin(sqrt(Pi)) # arcsine conversion
delta=4*var(theta.i)*(k-1) # sum of squares
N=2/(1/n1+1/n2) # harmonic mean, subcontrary mean
ncp1=N*delta # non-central parameter
power.chi=pchisq(qchisq(1-alpha,df),df,ncp1,lower=FALSE)
q.alpha=qchisq(1-alpha,df) # 3.841
qcrit=qchisq(max(p.val,1-p.val),df,ncp=0)
# qcrit = prop.test(c(r1,r2),c(n1,n2))$statistic
curve(dchisq(x,df),0,2*qcrit,xlab=quote(chi),ylab='Density',bty='n',lwd=2) # H0
curve(dchisq(x,df,ncp1),add=TRUE,lty=2,col=2,lwd.2) # H1
# power.chi=AUC of right half of the H1 curve
abline(v=qcrit)
abline(v=q.alpha,col='gray',lty=3)
legend('topright',bty='n',legend=c('H0','H1','chisq@p.value','chisq@alpha'),col=c(1,2,1,'gray'),lty=c(1,2,1,3),lwd=c(2,2,1,1))
text(qcrit,0,round(qcrit,2))
y0=dchisq(qcrit,df,0)
y1=dchisq(qcrit,df,ncp=ncp1)
FPR.equal=y0/(y0+y1) # length ratio
FPR.less=p.val/(p.val+power.chi) # area ratio
FPR.alpha=alpha/(alpha+power.chi) # FPR before analysis
VAL=c(power=power.chi,p.value=p.val,FPR.equal=FPR.equal,
FPR.less=FPR.less,FPR.alpha=FPR.alpha)
print(VAL,digits=3)
invisible(VAL)
}

499卵の名無しさん2018/05/19(土) 07:58:23.94ID:rrC4yXIM
PPV2prevalence <- function(sensitivity,specificity,PPV) {
(1-specificity)*PPV/((1-specificity)*PPV)+(1-PPV)*sensitivity))
}
FPP2prior <- function(power,FPR=0.05,pval=0.05){
pval*(1-FPR)/(pval*(1-FPR)+FPR*power)
}

500卵の名無しさん2018/05/19(土) 10:40:37.24ID:oTDRH91u
PPV2prevalence <- function(sensitivity,specificity,PPV) {
(1-specificity)*PPV/((1-specificity)*PPV+(1-PPV)*sensitivity)
}
PPV2prevalence(0.75,0.99,0.9)
FPP2prior <- function(power,FPR=0.05,pval=0.05){
pval*(1-FPR)/(pval*(1-FPR)+FPR*power)
}

501卵の名無しさん2018/05/20(日) 11:45:44.17ID:n2fbjQMc
# These date were used in 1908 by W. S. Gosset ('Student')
# as an example to illustrate the use of his t test,
# in the paper in which the test was introduced.
A=c(0.7,-1.6,-0.2,-1.2,-0.1,3.4,3.7,0.8,0.0,2.0)
B=c(1.9,0.8,1.1,0.1,-0.1,4.4,5.5,1.6,4.6,3.4)
t.test(A,B,var.equal = TRUE)
mean(A) ; mean(B)
sd(A) ; sd(B)
(E=mean(B)-mean(A))
nA=length(A) ; nB=length(B)
# The pooled estimate of the error within groups
SEpooled=sqrt(weighted.mean(c(var(A),var(B)),c(nA-1,nB-1)))
# Standard deviation of effect size
SE=sqrt(SEpooled^2/nA+SEpooled^2/nB)
(t.statistic=E/SE)
2*pt(t.statistic,nA+nB-2,lower=FALSE) # two-sided

502卵の名無しさん2018/05/20(日) 18:46:26.99ID:n2fbjQMc
# https://www.nejm.org/doi/10.1056/NEJMoa1207541
r1=14
n1=840
r2=73 # placebo
n2=829
calc.FPR.p2(r1,r2,n1,n2)
calc.FPR0.p2(r1,r2,n1,n2)
prop.test(c(r1,r2),c(n1,n2))


prior.needed <- function(r1,r2,n1,n2,FPR=0.05){
pval=prop.test(c(r1,r2),c(n1,n2))$p.value
ES=pwr::ES.h(r1/n1,r2/n2)
power=pwr::pwr.2p2n.test(ES,n1,n2,sig.level=pval)$power
prior=pval*(1-FPR)/(pval*(1-FPR)+FPR*power)
return(prior)
}

prior.needed(r1,r2,n1,n2,FPR=0.05)

503卵の名無しさん2018/05/21(月) 00:13:06.30ID:r/Pmsrg8
# https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4255087/#sup1
# p.508
# When p is a p-value with n1 samples, 95% ci of p-value of next experiment with n2 samples
# is supposed to be estimated by p2ci-function below.
p2ci <- function(p,n1=100,n2=100,sig.level=0.95){
lwr=pnorm(qnorm(p)*sqrt(n2/n1)-qnorm(1-(1-sig.level)/2)*sqrt(1+n2/n1))
# pnorm(qnorm(p1)-2.771808) # when n1=n2, 1.96 *s qrt(2) = 2.77
upr=pnorm(qnorm(p)*sqrt(n2/n1)+qnorm(1-(1-sig.level)/2)*sqrt(1+n2/n1))
# pnorm(qnorm(p1)+2.771808)
c(lwr,upr)
}
p2ci(0.05)
p2ci(0.001)
graphics.off()
pp=seq(0,1,length=1001)
plot(pp,sapply(pp,function(p)p2ci(p)[1]),type='l',bty='n',
xlab='initial p-value',ylab='C.I. of next p-value')
lines(pp,sapply(pp,function(p)p2ci(p)[2]))

###
(opt=optimise(function(p)p2ci(p)[1],c(0,1)))
p2ci(opt$minimum)
f <- function(n1=10,n2=10)c(t.test(rnorm(n1))$p.value,t.test(rnorm(n1))$p.value)
re=replicate(10^4,f())
points(re[1,],re[2,],col=rgb(0.01,0.01,0.01,0.01))

504卵の名無しさん2018/05/22(火) 15:20:40.05ID:pgTq+n72
等分散を仮定しないWelch法でのt検定からでも算出できるように

スクリプトを変更。

# https://www.physiology.org/doi/pdf/10.1152/advan.00042.2016
# Explorations in statistics: statistical facets of reproducibility Douglas Curran-Everett


p2p2w <- function(p.value,n1=10,n2=10,sd1=1,sd2=1,alpha=0.05){
p=p.value/2 # two.sided comparison
# Z-test
z=qnorm(p)
d.z=z*sqrt(sd1^2/n1+sd2^2/n2) # estimated difference of means by z.test
p2.z=pnorm(-qnorm(alpha/2)+z,lower=FALSE)
# Student
df=n1-1+n2-1
t=qt(p,df)
d.t=t*sqrt((1/n1+1/n2)*((n1-1)*sd1^2+(n2-1)*sd2^2)/((n1-1)+(n2-1)))
p2.t=pt(-qt(alpha/2,df)+t,df,lower=FALSE)
# Welch
var1=sd1^2 ; var2=sd2^2
dfw=(var1/n1+var2/n2)^2 / (var1^2/n1^2/(n1-1)+var2^2/n2^2/(n2-1))
t.w=qt(p,dfw)
d.w=t.w*sqrt(var1/n1+var2/n2)
p2.w=pt(-qt(alpha/2,dfw)+t.w,dfw,lower=FALSE)
data.frame(p2.w, p2.t, p2.z)
}

505卵の名無しさん2018/05/24(木) 21:03:43.35ID:R1FtDNpg
# simulation on the assumption of the event ~ poisson distribution
library(rstan)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)

stanString='
data{
int<lower=0> r1;
int<lower=0> r2;
int<lower=0> n1;
int<lower=0> n2;
}
parameters{
real<lower=0> lambda1;
real<lower=0> lambda2;
}
model{
r1 ~ poisson(lambda1) T[0,n1];
r2 ~ poisson(lambda2) T[0,n2];
}

'
# execute only on first occation
# stanModel=stan_model(model_code = stanString)
# saveRDS(stanModel,'2x2p.rds')

506卵の名無しさん2018/05/24(木) 21:05:13.78ID:R1FtDNpg
prop.poisson <- function(r1,r2,n1,n2, alpha=0.05){
data <- list(r1=r1,r2=r2,n1=n1,n2=n2)
stanModel=readRDS('2x2p.rds')
fit=sampling(stanModel,data=data,chain=1,iter=10000)
print(fit)
ms=rstan::extract(fit)
k=length(ms$lambda1)
p.sim=numeric(k)
for(i in 1:k){
p.sim[i]=prop.test(c(ms$lambda1[i],ms$lambda2[i]),c(n1,n2))$p.value
}
BEST::plotPost(p.sim)
mean(p.sim<alpha)
cat('Pr(p.value < ',alpha,') = ', mean(p.sim<alpha),'\n')
cat('Original p.value = ',prop.test(c(r1,r2),c(n1,n2))$p.value,'\n')
}

507卵の名無しさん2018/05/25(金) 22:02:04.22ID:8Ot6XASs

508卵の名無しさん2018/05/30(水) 10:36:17.99ID:5L7BTTj3
r=3
f<- function(x) (1-x)^(r-1)*x
curve (f(x))
auc=integrate (f,0,1)$value
pdf <- function (x) f(x)/auc
integrate(pdf, 0.5,1)$value
integrate (function (x)x*pdf(x),0,1)$value
2/(r+2)


z=3;N=9
f <- function (x) choose(N,z)*x^z*(1-x)^(N-z)
curve (f(x))
auc=integrate (f,0,1)$value
pdf <- function (x) f(x)/auc
integrate(pdf, 0.5,1)$value
integrate (function (x)x*pdf(x),0,1)$value
(z+1)/(N+2)

509卵の名無しさん2018/05/30(水) 20:54:56.95ID:w7KIivv+
f <- function(x) 1/sqrt(x*(1-x))
z <- function(x) integrate(f,0,x)$value
z(0.1)
z(0.5)
curve(f(x))
xx=seq(0,1,0.01)
plot(xx,sapply(xx,z))

510卵の名無しさん2018/05/30(水) 21:02:47.29ID:w7KIivv+
f <- function(x) 1/sqrt(x*(1-x))
z <- function(x) integrate(f,0,x)$value/pi
z(0.1)
z(0.5)
curve(f(x))
x=seq(.001,.999,by=.001)
plot(x,sapply(x,z))
abline(a=0,b=1)

511卵の名無しさん2018/05/30(水) 21:11:44.87ID:w7KIivv+
curve(pbeta(x,0.5,0.5))

512卵の名無しさん2018/05/30(水) 21:31:43.66ID:w7KIivv+
f <- function(x)dbeta(x,0.5,0.5)
zz <- function(x) integrate(f,0,x)$value
x=seq(.001,.999,by=.001)
plot(x,sapply(x,zz))
abline(a=0,b=1)

513卵の名無しさん2018/06/01(金) 11:46:08.93ID:UJHj2xQL
# marginal likelihoo of Fixed Effect Model
m.FE <- function(r1,n1,r2,n2,shape1=1,shape2=1){
choose(n1,r1)*choose(n2,r2)*beta(r1+r2+shape1,n1-r1+n2-r2+shape2)/beta(shape1,shape2)
}

# marginal likelihoo of Random Effect Model θ1~B(a1,ba),θ2~B(a2,b2)
m.REJ <- function(r1,n1,r2,n2,a1=1,b1=1,a2=1,b2=1){
2*choose(n1,r1)*choose(n2,r2)*integrate(
function(x){
dbeta(x,r1+a1,n1-r1+b1)*beta(r1+a1,n1-r1+b1)/beta(a1,b1)*
pbeta(x,r2+a2,n2-r2+b2,lower=FALSE)*beta(r2+a2,n2-r2+b2)/beta(a2,b2)
}
,0,1)$value
}

514卵の名無しさん2018/06/02(土) 14:59:28.21ID:iq0fCXui
Haldane <- function(x) 1/(x*(1-x))
curve(Haldane(x),0,1)
f <- function(x) integrate(Haldane,0.5,x)$value
x=seq(0.001,0.999,by=0.001)
plot(x,sapply(x,f),col=3)
curve(log(x/(1-x)),add=T)

515卵の名無しさん2018/06/04(月) 15:43:10.18ID:UVDFwwqx
エビデンス=周辺尤度= p(D |M1)=

p(D | θ,M1)p(θ |M1)dθ

516卵の名無しさん2018/06/04(月) 15:43:29.56ID:UVDFwwqx
エビデンス=周辺尤度
= p(D |M1)=

p(D | θ,M1)p(θ |M1)dθ

517卵の名無しさん2018/06/04(月) 18:46:43.68ID:UVDFwwqx
marginal likelihood of your model Mx is given by
p(D |Mx)= p(ξ1)p(D | ξ1)+p(ξ2)p(D | ξ2)+p(ξ3)p(D | ξ3)
=0.6×0.001+0.3×0.002+0.1×0.003
=0.0015.
The marginal likelihood is computed

518卵の名無しさん2018/06/04(月) 21:45:19.03ID:UVDFwwqx
stanStringJ='
data{
int<lower=0> N1;
int<lower=0> N2;
real Y1[N1];
real Y2[N2];
}
parameters{
real mu;
real<lower=0> sigma;
real alpha;
}
transformed parameters{
real mu1;
real mu2;
real delta;
real<lower=0> precision;
delta = alpha/sigma;
mu1 = mu + alpha/2;
mu2 = mu - alpha/2;
precision = 1/(sigma^2);
}
model{
Y1 ~ normal(mu1,sigma);
Y2 ~ normal(mu2,sigma);
mu ~ uniform(-9999,9999);
precision ~ gamma(0.5,0.5);
delta ~ cauchy(0,1);
}

'

519卵の名無しさん2018/06/04(月) 21:46:02.48ID:UVDFwwqx
N1=n1=137 ; m1=23.8 ; sd1=10.4 ; N2=n2=278 ; m2=25.2 ; sd2=12.3
U=scale(rnorm(n1))*sd1+m1 ; E=scale(rnorm(n2))*sd2+m2
Y1=as.vector(U) ; Y2=as.vector(E)
data=list(N1=N1,N2=N2,Y1=Y1,Y2=Y2)

# execute only on first occation
# JZS.model=stan_model(model_code = stanStringJ)
# JZS.model=stan_model('JZS.stan')
# saveRDS(JZS.model,'JZS.rds')
JZS.model=readRDS('JZS.rds')
fitJ=sampling(JZS.model,data=data, iter=20000)
print(fitJ,probs = c(.025,.5,.975))
ms=rstan::extract(fitJ)
dES=density(ms$delta)
plot(dES,lwd=2, xlab='Effect Size (δ)',main='The Savage?Dickey method',bty='l')
curve(dcauchy(x),lty=3, add=TRUE)
abline(v=0,col=8)
(d1=dES$y[which.min(dES$x^2)]) # posterior density at ES=0
(d0=dcauchy(0)) # prior density at ES=0
points(0,d0,cex=2)
points(0,d1,pch=19,cex=2)
legend('topleft',bty='n',legend=c('prior','posterior'),lty=c(3,1),lwd=c(1,2))
(BF10=d1/d0)
text(0,0,paste('BF10=',round(BF10,2)))

library(BayesFactor)
1/exp(ttestBF(Y1,Y2,rscale = 1)@bayesFactor[['bf']])

520卵の名無しさん2018/06/05(火) 09:10:30.95ID:Ob/ggqKh
N=10000
x=rbeta(N,0.5,0.5)
y=rbeta(N,0.5,0.5)
z=x-y
hist(z,freq=F)
dz=density(z)
dz$y[which.min(dz$x^2)]

z=x/y
hist(log10(z),freq=F)
dz=density(z)
min=1
dz$y[which.min((dz$x-min)^2)]

521卵の名無しさん2018/06/05(火) 09:14:46.94ID:Ob/ggqKh
data{ //binomBF.stan
int r1;
int r2;
int n1;
int n2;
real shape1;
real shape2;
}
parameters{
real <lower=0, upper=1> theta1;
real <lower=0, upper=1> theta2;
real <lower=0, upper=1> th1;
real <lower=0, upper=1> th2;
}
transformed parameters{
real rd;
real rr;
real rd0;
real rr0;
rd = theta1 - theta2;
rd0 = th1 - th2;
rr = theta1/theta2;
rr0 = th1/th2
}
model{
r1 ~ binomial(n1,theta1);
r2 ~ binomial(n2,theta2);
theta1 ~ beta(shape1,shape2);
theta2 ~ beta(shape1,shape2);
th1 ~ beta(shape1,shape2);
th2 ~ beta(shape1,shape2);
}

522卵の名無しさん2018/06/06(水) 21:19:25.17ID:IpiYYmjt
ある大学の入学者男女の比率は1であるという帰無仮説を検定する課題が花子と太郎に課された。

花子は50人を調査できたら終了として入学者を50人をみつけて18人が女子であるという結果を得た。
帰無仮説のもとで
50人中18人が女子である確率は 0.01603475
これ以下になるのは50人中0〜18人と32〜50人が女子の場合なので
両側検定して
> sum(dbinom(c(0:18,32:50),50,0.5))
[1] 0.06490865
> binom.test(18,50,0.5)$p.value
[1] 0.06490865
で帰無仮説は棄却できないと結論した。
http://i.imgur.com/XDIp9rM.png

523卵の名無しさん2018/06/06(水) 21:19:31.46ID:IpiYYmjt
一方、十八という数字が好きな太郎は一人ずつ調べて18人めの女子がみつかったところで調査を終えることにした。
18人めがみつかったのは花子と同じく50人めであった。
帰無仮説のもとで
18人がみつかるのが50人めである確率は0.005772512
これ以下になるのは23人以下50人以上番めで女子18人めがみつかった場合なので
両側検定して
pnb=dnbinom(0:999,18,0.5)
> 1 - sum(pnb[-which(pnb<=dnbinom(50-18,18,0.5))]) # < 0.05
[1] 0.02750309
http://i.imgur.com/K3T7utr.png
で帰無仮説は棄却される。

どちらの検定が正しいか、どちらも正しくないか?
検定する意図によってp値が変わるのは頻度主義統計の欠陥といえるか?

花子の横軸は女子数、太郎の横軸はサンプル数なので
サンプルでの女子の割合を横軸にして95%信頼区間を示す。
花子の検定での信頼区間は0.36〜0.72で18/50を含む、p=0.06491
http://i.imgur.com/SeTLk8K.jpg
太郎の検定での信頼区間は0.375〜0.72で18/50を含まない、p= 0.0275
http://i.imgur.com/tNzlfxe.jpg
主観である、検定の中止の基準の差でp値や信頼区間が変化するのは変だという批判である。

524卵の名無しさん2018/06/07(木) 07:36:58.84ID:6l5aJg03
50C18 * 0.5^18 * 0.5^32 と
49C17 * 0.5^17 * 0.5^32 * 0.5 の違いでしょう

18人目を見つけた人数を調べるというデザインがおかしいよね

これ事前確率0.5で50人調査して女が18人っていうのを
ベイズ更新していったらどうなる?

525卵の名無しさん2018/06/07(木) 07:47:39.40ID:Ni5wt/sw
>>524
酷いのになるとp<0.05になったらやめるとかいうのもあるな。
p-hackingと呼ばれる

526卵の名無しさん2018/06/08(金) 22:06:59.94ID:dTNUKiNw
r=21
N=20
a=0.5
b=0.5
p=a/(a+b+r+(1:N))
q=cumsum(p)
q
plot(1:N,p,ann=F)
plot(1:N,q,ann=F)


&#8203;

527卵の名無しさん2018/06/12(火) 00:02:23.51ID:XLL1LdWn
N=50
z=40
FP=0.01
shape1=1
shape2=1
data = list(N=N,z=z,FP=FP,shape1=shape1,shape2=shape2)
stanString=paste0('
data{
int N;
int z;
real FP;
real shape1;
real shape2;
}
parameters{
real<lower=0,upper=1> TP;
}
transformed parameters{
real<lower=0,upper=1> theta;
theta=((z*1.0)/(N*1.0)-FP)/(TP-FP);
}
model{
z ~ binomial(N,theta);
TP ~ beta(shape1,shape2);
}

')

528卵の名無しさん2018/06/12(火) 00:34:11.64ID:XLL1LdWn
N=50
z=40
FP=0.01
shape1=1
shape2=1
data = list(N=N,z=z,FP=FP,shape1=shape1,shape2=shape2)
stanString=paste0('
data{
int N;
int z;
real FP;
real shape1;
real shape2;
}
parameters{
real<lower=0,upper=1> TP;
}
transformed parameters{
real<lower=0,upper=1> theta;
theta=((z*1.0)/(N*1.0)-FP)/(TP-FP);
}
model{
z ~ binomial(N,theta);
TP ~ beta(shape1,shape2) T[0.5,];
}

')

model=stan_model(model_code = stanString)
fit=sampling(model,data=data,control=list(adapt_delta=0.99),iter=10000)
print(fit,digits=3,probs=c(.025,.50,.975))

529卵の名無しさん2018/06/12(火) 19:09:54.64ID:XLL1LdWn
# pdfからcdfの逆関数を作ってHDIを表示させて逆関数を返す
pdf2hdi <- function(pdf,xMIN=0,xMAX=1,cred=0.95,Print=TRUE){
nxx=1001
xx=seq(xMIN,xMAX,length=nxx)
xx=xx[-nxx]
xx=xx[-1]
xmin=xx[1]
xmax=xx[nxx-2]
AUC=integrate(pdf,xmin,xmax)$value
PDF=function(x)pdf(x)/AUC
cdf <- function(x) integrate(PDF,xmin,x)$value
ICDF <- function(x) uniroot(function(y) cdf(y)-x,c(xmin,xmax))$root
hdi=HDInterval::hdi(ICDF,credMass=cred)
print(c(hdi[1],hdi[2]),digits=5)
if(Print){
par(mfrow=c(3,1))
plot(xx,sapply(xx,PDF),main='pdf',type='h',xlab='x',ylab='Density',col='lightgreen')
legend('top',bty='n',legend=paste('HDI:',round(hdi,3)))
plot(xx,sapply(xx,cdf),main='cdf',type='h',xlab='x',ylab='Probability',col='lightblue')
pp=seq(0,1,length=nxx)
pp=pp[-nxx]
pp=pp[-1]
plot(pp,sapply(pp,ICDF),type='l',xlab='p',ylab='x',main='ICDF')
par(mfrow=c(1,1))
}
invisible(ICDF)
}

530卵の名無しさん2018/06/12(火) 19:12:46.61ID:XLL1LdWn
.n=50
.r=40
b=0.01
f = function(x,a,b,n,r){ # x:prevalence, a:TP, b:FP=0.01
p=x*a+(1-x)*b
choose(n,r)*p^r*(1-p)^(n-r)
}
f2 = function(x,n=.n,r=.r){
cubature::adaptIntegrate(function(ab)f(x,ab[1],b,n,r),
c(0,0),c(1,1))$integral
}
vf2=Vectorize(function(x)f2(x,n=.n,r=.r))
curve(vf2(x),ylab='',yaxs='i',axes=FALSE,lwd=2,xlab='prevalence') ; axis(1)
# points(0:10/10,vf2(0:10/10),pch=19)
optimise(function(x) vf2(x),c(0,1),maximum = TRUE)
auc=integrate(function(x)vf2(x),0,1)$value
pdf<-function(x)vf2(x)/auc
vpdf=Vectorize(pdf)
integrate(function(x)x*pdf(x),0,1)$value # mean
cdf=function(x)integrate(pdf,0,x)$value
vcdf=Vectorize(cdf)
# time consuming processes
# curve(vcdf(x),bty='l',n=64)
inv_cdf <- function(u){
uniroot(function(x,u0=u)vcdf(x)-u0,c(0,1),tol = 1e-18)$root
}
vinv_cdf=Vectorize(inv_cdf)
# curve(vinv_cdf(x),lwd=2,n=64)
hdi=HDInterval::hdi(vinv_cdf) ; hdi
# lower upper
# 0.7456652 1.0000000

531卵の名無しさん2018/06/12(火) 19:23:42.87ID:XLL1LdWn
# random numbers following PDF by John von Neuman's method
vonNeumann2 <- function(PDF,xmin=0,xmax=1,N=10000,Print=TRUE,...){
xx=seq(xmin,xmax,length=N+1)
xx=xx[-(N+1)]
xx=xx[-1]
ymax=max(PDF(xx))
Ux=runif(N,xmin,xmax)
Uy=runif(N,0,ymax)
Rand=Ux[which(Uy<=PDF(Ux))]
if(Print){
hist(Rand,xlim=c(xmin,xmax),freq=FALSE,col=sample(colors(),1),main='',...)
AUC=integrate(PDF,xmin,xmax)$value
lines(xx,sapply(xx,function(x)PDF(x)/AUC))
}
hdi=HDInterval::hdi(Rand)
print(c(hdi[1],hdi[2]),digits=5)
invisible(Rand)
}

532卵の名無しさん2018/06/12(火) 20:53:05.85ID:Ex3k8fq/
>>531
ここでもうりゅう先輩が迷惑掛けてんのか?

ウリュウなあ
こいつはなあ、生まれついてのビッグマウスであちこちに自分を売り込むが、
卒業しても国試浪人で医師免許ない50過ぎでは相手にされない
国試対策塾で非常識講師で細々と食つなぐが学生に馬鹿にされる
自分の医師コンプを隠すために医学生たちを「底辺」などという
実は自分が凄まじい底辺なのだが気づいていない

こんな嘘つきデブがのさばっているスレだな

ご苦労なこったよ、うりゅうのおっさん

わからねえとでも思ってんだろどーせ

533卵の名無しさん2018/06/13(水) 08:53:38.73ID:rcG4xYvl
library(rjags)
N=50
z=40
FP=0.01
shape1=1
shape2=1
dataList=list(N=N,z=z,FP=FP,shape1=shape1,shape2=shape2)
modelstring <- paste0("
model{
theta=TP*x+FP*(1-x)
z ~ dbinom(theta,N)
TP ~ dbeta(shape1,shape2)
x ~ dbeta(shape1,shape2)
}"
)
writeLines( modelstring , con="TEMPmodel.txt" )
jagsModel = jags.model( file="TEMPmodel.txt" , data=dataList, quiet=TRUE)
update(jagsModel)
codaSamples = coda.samples( jagsModel ,
variable=c("TP","x","theta"), n.iter=100000 )
js=as.matrix(codaSamples)
head(js)
BEST::plotPost(js[,'TP'],xlab='sensitivity')
BEST::plotPost(js[,'x'],xlab='prevalence')
BEST::plotPost(js[,'theta'],xlab='positive result',showMode = TRUE)

534卵の名無しさん2018/06/16(土) 16:23:17.07ID:V7A3qxKV
経理課の須藤は着服をやめろ!

勤務実態もないのに、グループ病院内から管理手当て(10万円)をもらうな!!!

意図的な給与操作、どうにかしろ!

535卵の名無しさん2018/06/17(日) 05:51:35.25ID:ifh2AARM
seqN <- function(N=100,K=5){
a=numeric(N)
for(i in 1:K) a[i]=2^(i-1)
for(i in K:(N-1)){
a[i+1]=0
for(j in 0:(K-1)){
a[i+1]=a[i+1]+a[i-j] # recursion formula
}
}

P0=numeric(N)
for(i in 1:N) P0[i]=a[i]/2^i # P0(n)=a(n)/2^n
P0

MP=matrix(rep(NA,N*K),ncol=K)
colnames(MP)=paste0('P',0:(K-1))
MP[,1]=P0
head(MP);tail(MP)
MP[1,2]=1/2
for(i in (K-2):K) MP[1,i]=0

for(k in 2:K){
for(i in 1:(N-1)) MP[i+1,k]=1/2*MP[i,k-1]
} # Pk(n+1)=1/2*P(k-1)(n)
ret=1-apply(MP,1,sum)

ret[N]
}
seqN(100,5)
seqN(1000,10)

536卵の名無しさん2018/06/17(日) 07:58:05.56ID:ifh2AARM
## 表の出る確率がpであるとき、N回コインを投げて K回以上表が連続する確率
seqNp <- function(N=100,K=5,p=0.5){
q=1-p
a=numeric(N) # a(n)=P0(n)/p^n , P0(n)=a(n)*p^n
for(i in 1:K) a[i]=q/p^i # P0(i)=q

for(i in K:(N-1)){ # recursive formula
a[i+1]=0
for(j in 0:(K-1)){
a[i+1]=(a[i+1]+a[i-j])
}
a[i+1]=q/p*a[i+1]
}

P0=numeric(N)
for(i in 1:N) P0[i]=a[i]*p^i # P0(n)=a(n)*p^n
P0

MP=matrix(rep(NA,N*K),ncol=K)
colnames(MP)=paste0('P',0:(K-1))
MP[,'P0']=P0
head(MP);tail(MP)
MP[1,'P1']=p
for(i in (K-2):K) MP[1,i]=0

for(k in 2:K){
for(i in 1:(N-1)) MP[i+1,k]=p*MP[i,k-1]
} # Pk(n+1)=p*P(k-1)(n)
ret=1-apply(MP,1,sum)

ret[N]
}

537卵の名無しさん2018/06/17(日) 08:33:04.03ID:LIEnVEKd
p:表
q=1-p
# Pk(n) (k=0,1,2,3,4)を途中、5連続して表が出ていなくて
# 最後のk回は連続して表が出ている確率とする。
#
P0(1)=q
P1(1)=p
P2(1)=P3(1)=P4(1)=0
P(k+1)(n+1)=p*Pk(n)
P0(n+1)=q*{P0(n)+P1(n)+P2(n)+P3(n)+P4(n)}
=q*{P0(n)+p*P0(n-1)+p^2*P0(n-2)+p^3*P0(n-3)+p^4*P0(n-4)}

P0(n)=a(n)*p^n
# a(n+1)p^(n+1)=q*p^n{a(n)+a(n-1)+a(n-2)+a(n-3)+a(n-4)}
# a(n+1)=q/p1*(a(n)+a(n-1)+a(n-2)+a(n-3)+a(n-4))

a(n)=P0(n)/p^n

538卵の名無しさん2018/06/17(日) 09:43:49.13ID:ifh2AARM
>>532

統計くらいできるのが国立卒の普通の臨床医。

おい、ド底辺

統計処理からはおまえは

都外のド底辺シリツ医大卒と推測されたが、あってるか?

539卵の名無しさん2018/06/17(日) 09:45:16.20ID:ifh2AARM
## 表の出る確率がpであるとき、N回コインを投げて K回以上表が連続する確率に一般化してみた。
seqNp <- function(N=100,K=5,p=0.5){
q=1-p
a=numeric(N) # a(n)=P0(n)/p^n , P0(n)=a(n)*p^n
for(i in 1:K) a[i]=q/p^i # P0(i)=q
for(i in K:(N-1)){ # recursive formula
a[i+1]=0
for(j in 0:(K-1)){
a[i+1]=(a[i+1]+a[i-j])
}
a[i+1]=q/p*a[i+1]
}

P0=numeric(N)
for(i in 1:N) P0[i]=a[i]*p^i # P0(n)=a(n)*p^n
MP=matrix(rep(NA,N*K),ncol=K)
colnames(MP)=paste0('P',0:(K-1))
MP[,'P0']=P0
head(MP);tail(MP)
MP[1,'P1']=p
for(i in (K-2):K) MP[1,i]=0

for(k in 2:K){
for(i in 1:(N-1)) MP[i+1,k]=p*MP[i,k-1]
} # Pk(n+1)=p*P(k-1)(n)
ret=1-apply(MP,1,sum)

ret[N]
}

540卵の名無しさん2018/06/17(日) 22:30:10.45ID:ifh2AARM
# pdfからcdfの逆関数を作ってHDIを表示させて逆関数を返す
pdf2hdi <- function(pdf,xMIN=0,xMAX=1,cred=0.95,Print=TRUE){
nxx=1001
xx=seq(xMIN,xMAX,length=nxx)
xx=xx[-nxx]
xx=xx[-1]
xmin=xx[1]
xmax=xx[nxx-2]
AUC=integrate(pdf,xmin,xmax)$value
PDF=function(x)pdf(x)/AUC
cdf <- function(x) integrate(PDF,xmin,x)$value
ICDF <- function(x) uniroot(function(y) cdf(y)-x,c(xmin,xmax))$root
hdi=HDInterval::hdi(ICDF,credMass=cred)
print(c(hdi[1],hdi[2]),digits=5)
if(Print){
par(mfrow=c(3,1))
plot(xx,sapply(xx,PDF),main='pdf',type='h',xlab='x',ylab='Density',col='lightgreen')
legend('top',bty='n',legend=paste('HDI:',round(hdi,3)))
plot(xx,sapply(xx,cdf),main='cdf',type='h',xlab='x',ylab='Probability',col='lightblue')
pp=seq(0,1,length=nxx)
pp=pp[-nxx]
pp=pp[-1]
plot(pp,sapply(pp,ICDF),type='l',xlab='p',ylab='x',main='ICDF')
par(mfrow=c(1,1))
}
invisible(ICDF)
}

541卵の名無しさん2018/06/18(月) 20:41:20.64ID:G5tabFnb
library(rjags)
modelstring <- paste0("
model{
theta=TP*x+FP*(1-x)
z ~ dbinom(theta,N)
TP ~ dbeta(shape1,shape2)
x ~ dbeta(shape1,shape2)
}"
)
writeLines( modelstring , con="TEMPmodel.txt" )
N=100 ; FP=0.01 ; shape1=1 ; shape2=1

guess.TP <- function(z){
dataList=list(N=N,z=z,FP=FP,shape1=shape1,shape2=shape2)
jagsModel = jags.model( file="TEMPmodel.txt" , data=dataList, quiet=TRUE)
update(jagsModel)
codaSamples = coda.samples(jagsModel,variable=c("TP","x"), n.iter=10000)
js=as.matrix(codaSamples)
m.TP=mean(js[,'TP'])
ci.TP=HPDinterval(as.mcmc(js[,'TP']))
m.x=mean(js[,'x'])
ci.x=HPDinterval(as.mcmc(js[,'x']))
c(m.TP=m.TP,ci.TP=ci.TP,m.x=m.x,ci.x=ci.x)
}
zz=1:20*5
re=sapply(zz,guess.TP)
head(re[,1:4])
re=as.matrix(re)
plot(zz,re['m.TP',],bty='l',ylim=c(0,1),type='n',las=1,
xlab='n : positives out of 100',ylab='sensitivity')
segments(zz,re[2,],zz,re[3,],col=8,lwd=3)
points(zz,re['m.TP',],pch=16)

542卵の名無しさん2018/06/18(月) 22:45:56.48ID:G5tabFnb
guess.TP2 <- function(z,FP){
dataList=list(N=N,z=z,FP=FP,shape1=shape1,shape2=shape2)
jagsModel = jags.model( file="TEMPmodel.txt" , data=dataList, quiet=TRUE)
update(jagsModel)
codaSamples = coda.samples(jagsModel,variable=c("TP"), n.iter=10000,thin=5)
js=as.matrix(codaSamples)
mean(js[,'TP'])
}
vG=Vectorize(guess.TP2)
n=1:20*5
FP=seq(0,0.25,length=20)
TP=outer(n,FP,vG) # wait several minutes
contour(n,FP,TP, col='navy',
xlab='n : positives out of 100',ylab='FP : 1-specificity',bty='l',nlevels=64)
points(50,0.10,pch='+',col='red',cex=1.5)

543卵の名無しさん2018/06/19(火) 22:41:00.16ID:ant1u+bV
n=5 # prime number
nn=1:(n-1)
tasu <- function(x,y) (x+y)%%n
hiku <- function(x,y) (x-y)%%n # row - col
kake <- function(x,y) (x*y)%%n
g=function(x) nn[which(x==1)]
.M=outer(nn,nn,kake)
G=apply(.M,2,g)
gyaku <- function(x) nn[which(G==(x%%n))]
waru <- function(x,y) (x*gyaku(y))%%n # row / col
waru(3,2)

xx=yy=c(0,nn)
names(xx)=paste0('x',c(0,nn))
names(yy)=paste0('y',c(0,nn))
outer(xx,yy,tasu) # x + y
outer(xx,yy,hiku) # x - y
outer(xx,yy,kake) # x * y

X=Y=nn
outer(X,Y,waru) # WRONG!!
outer(X,Y,Vectorize(waru))

a=expand.grid(X,Y)
b=matrix(mapply(waru,a[,1],a[,2]),ncol=length(X))
rownames(b)=paste0('x',nn)
colnames(b)=paste0('y',nn)
b # x / y

544卵の名無しさん2018/06/20(水) 19:25:59.38ID:myhyWcyK
rule3 <- function(n,confidence.level=0.95){
p=1/n
q=1-p # q^n.sample > 1-confidence.level
n.sample = log(1-confidence.level)/log(q)
return(n.sample)
}
n.sample=log(0.05)/log(0.999)

545卵の名無しさん2018/06/22(金) 14:47:07.52ID:ETsHYxXe
shuffle <- function(Cards){
n=length(Cards)
n1=ceiling(n/2)
n2=n-n1
C1=Cards[1:n1]
C2=Cards[(n1+1):n]
ret=NULL
for(i in 1:n1){
ret=c(ret,C1[i],C2[i])
}
ret[is.na(ret)==F]
}
x=as.character(c('A',2:10,'J','Q','K'))
cat(x,'\n') ; cat(shuffle(x))
Shuffles <- function(x){
tmp=shuffle(x)
i=1
while(!identical(x,tmp)){
tmp=shuffle(tmp)
i=i+1
}
return(i)
}

f =function(x)Shuffles(1:x)
nn=1:53
y=sapply(nn,f)
plot(nn,y,pch=16,bty='l',xlab='cards',ylab='shuffles')
cbind(nn,y)

546卵の名無しさん2018/06/22(金) 20:07:08.88ID:ETsHYxXe
http://000013.blogspot.com/2010/12/99.html

inversion <- function(x){ # 転倒数
n=length(x)
ret=numeric(n)
for(i in 1:(n-1)){
ret[i] = sum(x[i] > x[(i+1):n])
}
sum(ret)
}
x=c(4, 3, 5, 2, 1)
inversion(x)
is.even= function(x) !inversion(x)%%2
is.even(x)

prisoner99 <- function(n=100){
indx=sample(1:n,1) # defective number
X=sample((1:n)[-indx]) ; is.even(X)
Y=numeric(n-1)
for (i in 1:(n-1)){
x1=X[-i]
x2=(1:n)[!(1:n) %in% x1] # 囚人iが見えない番号
tmp=X
tmp[i]=x2[1] ; tmp[n]=x2[2]
Y[i]=ifelse(is.even(tmp), x2[1],x2[2]) # 偶順列になるように選択
}
all(X==Y)
}
mean(replicate(1e3,prisoner99()))

新着レスの表示
レスを投稿する