17/11/11 19:25:33.55 AQmX3Wf1.net
URLﾘﾝｸ(i.imgur.com)
251:卵の名無しさん
17/11/11 20:30:46.13 AQmX3Wf1.net
Golgo13 vs Golgo14
Let's MCMC with JAGS.
I set the ROPE ( Range of Practical Equivalence ) as -0.01 to 0.01
> smryMCMC( mcmcCoda2 , compVal=NULL ,ropeDiff = c(-0.01,0.01))
Mean Median Mode HDIlow HDIhigh CompVal PcntGtCompVal
Diff 0.07384483 0.05247579 0.01267166 -0.02478474 0.2345431 0 90.18
PcntLtROPE PcntInROPE PcntGtROPE
Diff 3.07 15.66 81.27
It'll be illustrated as follows.
URLﾘﾝｸ(i.imgur.com)
They are about equal to the result with Stan.
> mean(d)
[1] 0.07288324
> median(d)
[1] 0.05019718
> MAP(d)[1]# mode
0.01531325
> HDInterval::hdi(d)
lower upper
-0.0242909 0.2395352
> mean(d>0) # PcntGtCompVal
[1] 0.90575
> mean(d < -0.01) # PcntLtROPE
[1] 0.03475
> mean(-0.01 < d & d < 0.01) PcntInROPE
[1] 0.15975
> mean(d>0.01) # PcntGtROPE
[1] 0.8055
252:卵の名無しさん
17/11/11 23:28:55.02 AQmX3Wf1.net
0/5 vs 1/20
URLﾘﾝｸ(i.imgur.com)
253:卵の名無しさん
17/11/12 16:15:23.53 OKdyWAPj.net
posterior ∝ likelihood * prior
URLﾘﾝｸ(i.imgur.com)
254:卵の名無しさん
17/11/13 10:06:43.47 9LmqAQrY.net
K=12
omega1=0.25
omega2=0.75
a1 = omega1*(K-2)+1
b1 = (1-omega1)*(K-2)+1
a2 = omega2*(K-2)+1
b2 = (1-omega2)*(K-2)+1
curve(dbeta(x,a1,b1))
curve(dbeta(x,a2,b2))
curve(dbeta(x,a1,b1)+dbeta(x,a2,b2))
pdf <- function(theta) (dbeta(theta,a1,b1)+dbeta(theta,a2,b2))/2
n=31
d=0.02
h=function(theta,omega){ # some omega must equal to omega1 or omega2
if(omega>=omega1-d & omega<=omega1+d){
a1 = omega*(K-2)+1
b1 = (1-omega)*(K-2)+1
return(dbeta(theta,a1,b1))
}
if(omega>=omega2-d & omega<=omega2+d){
a2 = omega*(K-2)+1
b2 = (1-omega)*(K-2)+1
return(dbeta(theta,a2,b2))
}
return(0)
}
255:卵の名無しさん
17/11/13 10:07:02.25 9LmqAQrY.net
prior=matrix(NA,n,n) # with outer ERROR!
for(i in 1:n){
for(j in 1:n){
prior[i,j]=h(theta[i],omega[j])
}
}
image(theta,omega,prior,col = heat.colors(12,0.5))
contour(theta,omega,prior,add=TRUE)
persp(theta,omega,prior,col='skyblue',theta = 30)
library(rgl)
persp3d(theta,omega,prior,col='skyblue')
s=6 ; f=3
h2=function(theta,omega) theta^s*(1-theta)^f
likelihood=outer(theta,omega,h2)
image(theta,omega,likelihood,col = terrain.colors(12,0.5))
contour(theta,omega,likelihood,add=TRUE)
persp(theta,omega,likelihood,col='wheat',theta = 30)
library(rgl)
persp3d(theta,omega,likelihood,col='wheat')
posterior=prior*likelihood
image(theta,omega,posterior,col = topo.colors(12,0.5))
contour(theta,omega,posterior,add=TRUE)
persp(theta,omega,posterior,col='maroon',theta = 30)
library(rgl)
persp3d(theta,omega,posterior,col='maroon')
256:卵の名無しさん
17/11/13 12:14:38.18 9LmqAQrY.net
ド底辺シリツ医大は裏口入学と学力で入った例外入学がいるとする。
高卒レベルの基礎学力テストをしたところ裏口入学は不合格率の最頻値が0.75、例外者のそれは0.25であった。
いずれの分布も形状母数和が12のベータ分布に従っていた。
ド底辺シリツ医大でテストしたところ6人が不合格、3人が合格であったとき、ド底辺シリツ医大の裏口入学者の割合を推測せよ。
257:卵の名無しさん
17/11/13 20:42:17.05 JnL2ZVyT.net
omega1=0.25
omega2=0.75
K=12
modelString='
model {
f o r(ii n1 : N){
y[i] ~ dbern( theta )
}
theta ~ dbeta( omega[m]*(kappa-2)+1 , (1-omega[m])*(kappa-2)+1 )
omega[1] <-
omega1
omega[2] <- omega2
kappa <- K
m ~ dcat( mPriorProb[] )
mPriorProb[1] <- p
mPriorProb[2] <- 1 - p
p ~ dunif(0,1)
}
'
dataList=list(omega1=omega1,omega2=omega2,K=K)
writeLines( modelString , con="TEMPmodel.txt" )
jagsModel = jags.model( file="TEMPmodel.txt" , data=dataList)
update( jagsModel)
codaSamples = coda.samples( jagsModel , variable=c('omega','kappa','theta'), n.iter=10000)
summary(codaSamples)
plot(codaSamples)
258:卵の名無しさん
17/11/13 20:51:13.12 JnL2ZVyT.net
# in debug
y= c(1,1,1,1,1,1,0,0,0)
N=length(y)
omega1=0.25
omega2=0.75
K=12
modelString='
model {
for(i n1:N){
y[i] ~ dbern( theta )
}
theta ~ dbeta( omega[m]*(kappa-2)+1 , (1-omega[m])*(kappa-2)+1 )
omega[1] <- omega1
omega[2] <- omega2
kappa <- K
m ~ dcat(mPriorProb[] )
mPriorProb[1] <- p
mPriorProb[2] <- 1 - p
p ~ dunif(0,1)
}
'
dataList=list(y=y,N=N, omega1=omega1,omega2=omega2,K=K)
writeLines( modelString , con="TEMPmodel.txt" )
jagsModel = jags.model( file="TEMPmodel.txt" , data=dataList)
update( jagsModel)
codaSamples = coda.samples( jagsModel , variable=c('m,'theta'), n.iter=10000)
summary(codaSamples)
plot(codaSamples)
259:卵の名無しさん
17/11/13 21:18:42.89 JnL2ZVyT.net
# in debug
y= c(1,1,1,1,1,1,0,0,0)
N=length(y)
omega1=0.25
omega2=0.75
kappa1=12
kappa2=12
modelString='
model {
for (i in 1:N){
y[i] ~ dbern( theta )
}
theta <- equals(m,1)*theta1 + equals(m,2)*theta2
theta1 ~ dbeta( omega1*(kappa1-2)+1 , (1-omega1)*(kappa1-2)+1 )
theta2 ~ dbeta( omega2*(kappa2-2)+1 , (1-omega2)*(kappa2-2)+1 )
m ~ dcat( mPriorProb[] )
mPriorProb[1] <- p
mPriorProb[2] <- 1-p
p ~ dunif(0,1)
}
'
dataList=list(y=y,N=N, omega1=omega1,omega2=omega2, kappa1=kappa1, kappa2=kappa2)
writeLines( modelString , con="TEMPmodel.txt" )
jagsModel = jags.model( file="TEMPmodel.txt" , data=dataList)
update( jagsModel)
codaSamples = coda.samples( jagsModel , variable=c('m,'theta'), n.iter=10000)
summary(codaSamples)
plot(codaSamples)
260:卵の名無しさん
17/11/13 21:51:39.83 JnL2ZVyT.net
P(Uraguchi|Wrong) = P(Wrong|Uraguchi)*P(Uraguchi)/(P(Wrong|Uraguchi)*P(Uraguchi)+P(Wrong|Square)*P(Square))
Uraguchi: Buying the way
261:to Do-Teihen, unfair matriculation Wrong: Writing Wrong English Square: fair matriculation P(Square) = 1 - P(Uraguchi) P(Wrong|Uraguchi) = 1 P(Wrong|Square)=0.001 P(Uraguchi) ~ dunif(0,1) P(U|W)=1*P(U)/(1*P(U)+0.001*(1-P(U))) modelString=' puw <- 1 * pu /(1 * pu + 0.001*(1-pu)) pu ~ dunif(0,1) '
262:卵の名無しさん
17/11/14 07:16:26.27 kN15uX//.net
>>249
> js=as.matrix(codaSamples)
> boxplot(theta~m)
> tapply(theta,m,length)
1 2
8778 41222
> sum(m==2)/length(m)
[1] 0.82444
263:卵の名無しさん
17/11/14 11:23:12.95 /SJKjWLk.net
require(rjags)
JmodelString='
model {
puw = 1 * pu /(1 * pu + 0.001*(1-pu))
pu ~ dunif(0,1)
}
'
writeLines( JmodelString , con="TEMPmodel.txt" )
jagsModel = jags.model( file="TEMPmodel.txt")
update( jagsModel)
codaSamples = coda.samples( jagsModel , variable=c('pu','puw'), n.iter=100000)
summary(codaSamples)
js=as.matrix(codaSamples)
hist(js[,'puw'],xlim=c(0.9,1),freq=FALSE,breaks=1000,col='red',
main='P(Uraguchi|Wrong_English)',xlab='Probability')
HDInterval::hdi(js[,'puw'])
264:卵の名無しさん
17/11/14 13:59:12.73 /SJKjWLk.net
# Stan for Uraguchi
modelString='
parameters{
real<lower=0,upper=1> pu;
}
transformed parameters{
real puw = 1 * pu /(1 * pu + 0.001*(1-pu));
}
model{
pu ~ uniform(0,1);
}
ura.model<-stan_model(model_code = modelString)
fit.ura <- rstan::sampling(ura.model,seed=123,iter=10000,warmup=5000)
print(fit.ura,digits=4)
ms=rstan::extract(fit.ura)
round(HDInterval::hdi(ms$puw),4)
265:卵の名無しさん
17/11/14 18:44:34.87 Njwo3HI0.net
ド底辺シリツ医大は裏口入学と学力で入った例外入学がいるとする。
高卒レベルの基礎学力テストをしたところ裏口入学は不合格率の最頻値が0.75、例外者のそれは0.25であった。
いずれの分布も形状母数和が12のベータ分布に従っていた。
あるド底辺シリツ医大でテストしたところ6人が不合格、3人が合格であったとき、このド底辺シリツ医大の裏口入学者の割合を推測せよ。
Suppose there are Uraguchi and non-Uraguchi(irregular) enrollees in the DoTeihen medical school,
they get the achievement test for high school students.
The failing rate of Uraguchi is known to distribute in β distributuion with its mode value ω = 0.75 and sum of shape parameters κ = 12.
The failing rate of irregulars is in β distribution with ω = 0.25 and κ = 12.
At a DoTeihen medical school, 9 alimni had the test, and 6 failed and 3 passed.
Infer the proportion of Uraguchi in this DoTeihen.
266:卵の名無しさん
17/11/15 10:17:36.22 4qVYF7j+.net
y= c(1,1,1,1,1,1,1,1,1)
N=length(y)
omega1=0.75
omega2=0.25
kappa1=12
kappa2=12
modelString='
model {
for (i in 1:N){
y[i] ~ dbern( theta )
}
theta <- equals(m,1)*theta1 + equals(m,2)*theta2
theta1 ~ dbeta( omega1*(kappa1-2)+1 , (1-omega1)*(kappa1-2)+1 )
theta2 ~ dbeta( omega2*(kappa2-2)+1 , (1-omega2)*(kappa2-2)+1 )
m ~ dcat( mPriorProb[] )
mPriorProb[1] <- p
mPriorProb[2] <- 1-p
p ~ dunif(0,1)
}
'
dataList=list(y=y,N=N, omega1=omega1,omega2=omega2, kappa1=kappa1, kappa2=kappa2)
writeLines( modelString , con="TEMPmodel.txt" )
jagsModel = jags.model( file="TEMPmodel.txt" , data=dataList)
update( jagsModel)
codaSamples = coda.samples( jagsModel , variable=c('m','theta','theta1','theta2'), n.iter=50000)
summary(codaSamples)
js=as.matrix(codaSamples)
tapply(js[,'theta'],js[,'m'],length)
sum(js[,'m']==1)/nrow(js)
267:卵の名無しさん
17/11/15 14:20:17.72 RUmagzE5.net
In Bayesian data analysis, evidence is the marginal likelihood （Integrate P(D|Θ)P(Θ)dΘ） which MCMC cannot yield.
268:卵の名無しさん
17/11/17 12:56:37.51 nii6SWM6.net
# randam numbers by inverse culmutive density function
RandICDF <- function(ICDF,PDF,...){
U=runif(10000)
rand=ICDF(U,...)
hist(rand,freq=FALSE,breaks=30,
col=sample(colors(),2),main='')
curve(PDF(x,...),add=TRUE,lwd=2)
invisible(rand)
}
par(mfrow=c(2,2))
RandICDF(qnorm,dnorm)
RandICDF(qbeta,dbeta,shape1=3.5,shape2=8.5)
RandICDF(qbeta,dbeta,shape1=0.1,shape2=0.1)
RandICDF(qgamma,dgamma,shape=3,rate=1)
URLﾘﾝｸ(i.imgur.com)
269:卵の名無しさん
17/11/17 13:24:40.98 nii6SWM6.net
# randam numbers following PDF by Jhon von Neuman's method
vonNeumann <- function(PDF,xmin=0,xmax=1){
N=10000
ymax=max(PDF(seq(xmin,xmax,length=N+1)))
Ux=runif(N,xmin,xmax)
Uy=runif(N,0,ymax)
Rand=Ux[which(Uy<=PDF(Ux))]
hist(Rand,xlim=c(xmin,xmax),freq=FALSE,breaks=30,col=sample(colors(),2),main='')
curve(PDF,add=TRUE,lwd=2)
invisible(Rand)
}
par(mfrow=c(2,2))
vonNeumann(function(x)dbeta(x,3.5,8.5))
vonNeumann(function(x)dgamma(x,3,1),0,10)
vonNeumann(dnorm,xmin=-3,xmax=3)
p=runif(1)
f=function(x) p*dnorm(x,-3,1)+(1-p)*dnorm(x,3,3)
vonNeumann(f,xmin=-10,xmax=10)
URLﾘﾝｸ(i.imgur.com)
270:卵の名無しさん
17/11/17 16:06:01.69 nii6SWM6.net
Golgo Script will be reduced as follows:
N shoots with z hits
Nz <- function(N,z,...){
curve(dbeta(x,1+z,1+N-z),xlab='Probability of Hit',
ylab=paste('Probabilty of',z,'
271:Hits out of ',N,'Shoots'),...) hdi=HDInterval::hdi(qbeta,shape1=1+z,shape2=1+N-z) print(hdi,digits=4) }
272:卵の名無しさん
17/11/18 10:52:47.33 V/eIhsOZ.net
modelString =
'model {
z ~ dbin(0.5,N)
N ~ dpois(lambda)
p=z/N
}
'
writeLines(modelString , con="TEMPmodel.txt" )
f <- function (lambda){
dataList=list(lambda=lambda)
jagsModel = jags.model( file="TEMPmodel.txt" , data=dataList)
update( jagsModel)
codaSamples = coda.samples( jagsModel , variable=c('p'), n.iter=10000)
js=as.matrix(codaSamples)
jsp=js[,'p']
mean(jsp<=7/24)
}
xx=7:50
yy=sapply(xx,f)
plot(xx,yy)
273:卵の名無しさん
17/11/18 16:59:45.52 FiDbN6uZ.net
graphics.off()
p=binom.test(7,24,alt='less')$p.value
1-(1-p)^2 # 0.06289339
z1=7
N1=24
N2=24
plot(0:N1/N1,dbinom(0:N1,N1,0.5),type='h',lwd=5,col='skyblue',ann=FALSE,ylim=c(0,0.25))
points(z1/N1,0,pch='+',cex=2)
binom.test(z1,N1,alt='less')$p.value
0:N1/N1 < z1/N1
sum(0:N1/N1 <= z1/N1) # 8個0:7
sum(dbinom(0:z1,N1,0.5))
# A∪B == A + B - A∩B 0.06289339
sum(dbinom(0:7,N1,0.5))+sum(dbinom(0:7,N2,0.5)) -sum(outer(dbinom(0:7,N1,0.5),dbinom(0:7,N2,0.5),'*'))
z1=7
N1=24
N2=12
plot(0:N2/N2,dbinom(0:N2,N2,0.5),type='h',lwd=5,col='wheat',ann=FALSE)
lines(0:N1/N1+0.05,dbinom(0:N1,N1,0.5),type='h',lwd=5,col='skyblue',ann=FALSE,ylim=c(0,0.25))
points(z1/N1,0,pch='+',cex=2)
0:N2/N2 < z1/N1
sum(0:N2/N2 <= z1/N1) # 4個0:3
# A∪B == A + B - A∩B 0.1026226
sum(dbinom(0:7,N1,0.5))+sum(dbinom(0:3,N2,0.5)) - sum(outer(dbinom(0:z1,N1,0.5),dbinom(0:3,N2,0.5)))
274:卵の名無しさん
17/11/20 19:17:52.08 xJug4kDO.net
#
pdf2hdi <- function(pdf,xMIN=0,xMAX=1,cred=0.95,Print=FALSE){
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
cdf <- function(x) integrate(pdf,xmin,x)$value/AUC
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]))
if(Print){
pp=seq(0,1,length=nxx)
plot(pp,sapply(pp,ICDF),type='l',lwd=2,xlab='p',ylab='x')
legend('top',bty='n',legend=paste('HDI:',round(hdi,3)))
}
invisible(ICDF)
}
pdf2hdi(function(x)dbeta(x,0.001,0.001)*x^7*(1-x)^17)
pdf2hdi(dnorm,-10,10,cre=0.95)
275:卵の名無しさん
17/11/20 20:38:04.94 xJug4kDO.net
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=4)
if(Print){
par(mfrow=c(3,1))
plot(xx,sapply(xx,PDF),main='pdf',type='l',lwd=2,xlab='x',ylab='Density')
legend('top',bty='n',legend=paste('HDI:',round(hdi,3)))
plot(xx,sapply(xx,cdf),main='cdf',type='l',lwd=2,xlab='x',ylab='Probability')
pp=seq(0,1,length=nxx)
pp=pp[-nxx]
pp=pp[-1]
plot(pp,sapply(pp,ICDF),type='l',lwd=2,xlab='p',ylab='x',main='ICDF')
par(mfrow=c(1,1))
}
invisible(ICDF)
}
276:卵の名無しさん
17/11/20 21:23:09.45 xJug4kDO.net
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)
}
277:卵の名無しさん
17/11/20 22:09:25.75 xJug4kDO.net
According to this posting, ｽﾚﾘﾝｸ(doctor板:529番) as many as 15 freshmen flunk.
Let's assume there are 120 students in one class and the number of flunker is distributed as poisson distribution.
In what range will students are expected to flunk next year? Calculate the number with 99% confidence interval.
pdf2hdi <- function(pdf,xMIN=0,xMAX=1,cred=0.95,Print=TRUE,FUN=FALSE){
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))
}
if(FUN) invisible(ICDF)
invisible(hdi)
}
278:卵の名無しさん
17/11/22 10:27:31.38 WHgLvp/3.net
ab2mv<-function(a,b){
m<-a/(a+b)
v<-m*(1-m)/(a+b+1)
mv<-c(m,v)
return(mv)
}
mv2ab<-function(m,v){
a=(-m^3+m^2-m*v)/v
b=(m^3-2*m^2+m*v+m-v)/v
ab<-c(a,b)
return(ab)
}
HDCI <- function(PMF,cl=0.95){ # Highest Density Confidence Interval
PDF=PMF/sum(PMF)
rsPDF=rev(sort(PDF))
min.density=rsPDF[min(which(cumsum(rsPDF)>=cl))]
index=which(PDF>=min.density)
data.frame(lower.idx=round(min(index)),upper.idx=round(max(index)),actual.CI=sum(PDF[index]))
}
279:卵の名無しさん
17/11/22 12:17:28.03 WHgLvp/3.net
Scale <- function(x) (x-mean(x))/sqrt(sum((x-mean(x))^2)/(length(x)-1))
m=50;s=10
rn
280:=rnorm(1000,m,s) mean(rn) sd(rn) y=Scale(rn)*s+m mean(y) sd(y)
281:卵の名無しさん
17/11/22 19:57:10.44 WHgLvp/3.net
URLﾘﾝｸ(youtu.be)
282:卵の名無しさん
17/11/23 15:22:14.60 Ue5tZuwc.net
par(mfrow=c(2,2))
theta=0.5
NN=1000
N=1:NN
flip=numeric()
for(i in N){
flip=append(flip,rbinom(1,1,theta))
}
z=cumsum(flip)
z_N=z/N
plot(z_N,type='l',ylim=c(0,1))
pv=numeric()
for(i in N){
pv[i]=binom.test(z[i],N[i],theta)$p.value
}
plot(pv,type='l')
abline(h=0.05,col='blue',ylim=c(0,1))
bf=numeric()
for(i in N){
bf[i]=beta(z[i]+1,N[i]-z[i]+1)/beta(1,1)/(theta^z[i]*(1-theta)^(N[i]-z[i]))
}
plot(log(bf),type='l')
abline(h=log(3),lty=3)
abline(h=log(1/3),lty=3)
283:卵の名無しさん
17/11/23 15:22:43.97 Ue5tZuwc.net
hdi=NULL
for(i in N){
y=flip[1:i]
s=rep(1,i)
data=data.frame(y,s)
Ntotal=i
Nsubj=1
dataList=list(y=y,s=s,Ntotal=Ntotal,Nsubj=Nsubj)
js=genMCMC2(data)
hdi=rbind(hdi,HDIofMCMC(as.matrix(js)))
}
saveRDS(hdi,'hdi_sequential')
hdi=readRDS('hdi_sequential')
plot(hdi[,1],type='l',ylim=c(0,1),main='95% HDI')
lines(hdi[,2])
abline(h=0.5,col=4)
plot(apply(hdi,1,diff),type='l',main='HDI width')
284:卵の名無しさん
17/11/25 12:10:33.25 3kOACkIe.net
Metropolis Algo
source('DBDA2E-utilities.R')
Bern <- function(x) rbinom(1,1,x)
Metro <- function(
SD=0.02,
.z=14,
.N=20,
a=1,
b=1,
k=50000,
Print=TRUE){
theta=numeric()
theta[1]=0.01
likely <- function(x,z=.z,N=.N) x^z*(1-x)^(N-z)
for(i in 1:(k-1)){
delta=rnorm(1,0,SD)
theta_p=theta[i]+delta
p=min(1,likely(theta_p)*dbeta(theta_p,a,b)/
(likely(theta[i])*dbeta(theta[i],a,b)))
theta[i+1]=theta[i]+Bern(p)*delta
}
if(Print){
# plot(theta,1:k,type='l')
plotPost(theta,cenTend = 'mean',cex.lab = 1)
plot(theta[1:100],1:100,type='l',xlim=c(0,1))
plot(theta[(k-100):k],(k-100):k,type='l',xlim=c(0,1))}
invisible(theta)
}
285:卵の名無しさん
17/11/26 19:08:27.54 11LwtYgG.net
# JAGS for proportion
graphics.off()
rm(list=ls())
zi=100; Ni=100; zj=10; Nj=10
(y=c(rep(1,zi),rep(0,Ni-zi),rep(1,zj),rep(0,Nj-zj)))
(s=as.numeric(factor(c(rep('D',Ni),rep('U',Nj)))))
myData=data.frame(y=y,s=s)
Ntotal = length(y)
Nsubj = length(unique(s))
dataList = list(
y = y ,
s = s ,
Ntotal = Ntotal ,
Nsubj = Nsubj
)
a=1 ; b=1 # prior : beta(a,b)
# JAGS model
modelString = paste0("
model {
for ( i in 1:Ntotal ) {
y[i] ~ dbern( theta[s[i]] )
}
for ( sIdx in 1:Nsubj ) {
theta[sIdx] ~ dbeta(", a,',' , b," )
}
}
")
# close modelString
286:卵の名無しさん
17/11/26 19:09:02.05 11LwtYgG.net
writeLines( modelString , con="TEMPmodel.txt" )
jagsModel = jags.model( file="TEMPmodel.txt" , data=dataList)
#jags.model(file, data, inits,n.chains = 1, n.adapt=1000, quiet=FALSE)
update( jagsModel)
codaSamples = coda.samples( jagsModel , variable="theta", n.iter=10000 )
#coda.samples(model, variable.names, n.iter, thin = 1, na.rm=TRUE, ...)
summary(codaSamples)
plot(codaSamples,col=sample(colors()))
mcmcMat=as.matrix(codaSamples)
par(mfrow=c(2,2))
hist(mcmcMat[,1],freq=FALSE,xlim=c(0,1),
col=sample(colours(),1))
hist(mcmcMat[,1]-mcmcMat[,2],freq=FALSE,xlim=c(-1,1),
col=sample(colours(),1))
plot(mcmcMat[,1],mcmcMat[,2], col=rgb(0.01,0.01,0.3,0.25))
hist(mcmcMat[,2],freq=FALSE, xlim=c(0,1),
col=sample(colours(),1))
287:卵の名無しさん
17/11/26 19:09:08.75 11LwtYgG.net
source("Kruschke_tools.R") # for genMCMC2
(myData=data.frame(y,s))
mcmcCoda = genMCMC2( data=myData , numSavedSteps=10000,a=1,b=1)
# mcmcCoda = genMCMC( data=myData , numSavedSteps=10000)
# genMCMC
# Display diagnostics of chain, for specified parameter:
source('Jags-Ydich-XnomSsubj-MbernBeta.R')
diagMCMC( mcmcCoda , parName="theta[1]" )
diagMCMC( mcmcCoda , parName="theta[2]" )
# Display numerical summary statistics of chain:
smryMCMC( mcmcCoda , compVal=NULL,ropeDiff = c(-0.025,0.025))
#function( codaSamples , compVal=NULL , rope=NULL , saveName=NULL )
summary(mcmcCoda)
# Display graphical posterior information:
plotMCMC( mcmcCoda , data=myData , compVal=NULL,ropeDiff = c(-0.025,0.025))
# function( codaSamples , data , compVal=NULL , rope=NULL ,
# saveName=NULL , showCurve=FALSE , saveType="jpg" )
plot(mcmcCoda)
plotMCMC2( mcmcCoda , data=myData , compVal=NULL, showCurve=FALSE,
.credMass = 0.95,ropeDiff = c(-0.025,0.025),cenTend='mean')
(summry=smryMCMC( mcmcCoda, compVal=NULL, ropeDiff = c(-0.025,0.025)))
print(summry[3,],digits=2)
prop.test(c(zi,zj),c(Ni,Nj))$p.value
fisher.test(matrix(c(zi,Ni-zi,zj,Nj-zj),2))$p.value
288:卵の名無しさん
17/11/28 18:49:23.50 QFCizNPN.net
logistic <- function(x,gain,threshol
289:d){ 1/(1 + exp(-gain*(x-threshold))) } b2gt <- function(b0,b1,b2){ gain=sqrt(b1^2+b2^2) threshold=-b0/gain return(gain=gain,threshold=threshold) }
290:卵の名無しさん
17/11/28 18:52:07.55 QFCizNPN.net
I could not locate a good site to explain normalizationn for logististic regression,
but with the examples depicted in the textbook I have finally got understood.
This is it.
URLﾘﾝｸ(i.imgur.com)
core portion of its code :
ｽﾚﾘﾝｸ(hosp板:275番)
291:卵の名無しさん
17/11/28 20:13:41.22 QFCizNPN.net
# z=b0+巴k*xk
# p = logistic(z) = 1/(1+e^-z)
# 1-p = (1+e^-z)/(1+e^-z) - 1/(1+e^-z) = e^-z/(1+e^-z)
# p/(1-p) = 1 /e^-z = e^z
# logit(p) = log(p/(1-p))= log(e^z) = z
# logit(logistic(z)) = z
292:卵の名無しさん
17/11/30 16:34:37.18 TyAFrmPC.net
dataList=list(y=y,Ntotal=length(y),meanY=mean(y),sdY=sd(y))
modelString = '
model {
for ( i in 1:Ntotal ) {
y[i] ~ dt( mu , 1/sigma^2 , nu )
}
mu ~ dnorm( meanY , 1/(100*sdY)^2 )
sigma ~ dunif( sdY/1000 , sdY*1000 )
nu ~ dexp(1/30.0)
}
'
writeLines(modelString , con="TEMPmodel.txt" )
jagsModel = jags.model( file="TEMPmodel.txt" , data=dataList)
update( jagsModel)
codaSamples = coda.samples( jagsModel , variable=c('sigma','mu','nu'), n.iter=10000)
plot(codaSamples,col=sample(colours(),1))
js=as.matrix(codaSamples)
head(js)
293:卵の名無しさん
17/11/30 16:35:08.68 TyAFrmPC.net
# Y = aX + b , X ~ dt, a:scale parameter, b:location parameter
dt_ls <- function(x, df, mu, a) 1/a * dt((x - mu)/a, df)
pt_ls <- function(x, df, mu, a) pt((x - mu)/a, df)
qt_ls <- function(prob, df, mu, a) qt(prob, df)*a + mu
rt_ls <- function(n, df, mu, a) rt(n,df)*a + mu
par(mfrow=c(1,1))
hist(y,breaks=20,col='skyblue',freq=FALSE,xlim=c(30,220),main='')
N=63 #length(y)
for(i in sample(1:nrow(js),N)){
curve(dt_ls(x,js[i,'nu'],js[i,'mu'],js[i,'sigma']),add=TRUE,
lty=1,col=rgb(.01,.01,.01,.1))
}
294:卵の名無しさん
17/12/01 13:08:02.18 UfpWtEOZ.net
dT <- function(x, nu, mu, sd){
s=sd*sqrt((nu-2)/nu)
dt((x - mu)/s, nu)/s
}
pT <- function(x, nu, mu, sd){
s=sd*sqrt((nu-2)/nu)
pt((x - mu)/a, nu)
}
qT <- function(prob, nu, mu, sd){
s=sd*sqrt((nu-2)/nu)
qt(prob, nu)*s + mu
}
rT <- function(n, nu, mu, sd){
s=sd*sqrt((nu-2)/nu)
rt(n,nu)*s + mu
}
295:卵の名無しさん
17/12/02 06:06:16.06 SDqtqHE2.net
男性 28.2％
女性 9.0％
男女計 18.2％
URLﾘﾝｸ(www.jti.co.jp)
P(s)=.182
P(s|f)=.090
P(s|m)=.282
P(s)=P(s|f)P(f)+P(s|m)P(m)
P(m)=1-P(f)
から
P(f)=(P(s)-P(s|m))/(P(s|f)-P(s|m))
ベイズの公式
P(f|s)=P(s|f)P(f)/(P(s|m)P(m)+P(s|f)P(f))
P(s|m)P(m)+P(s|f)P(f)=P(s)
.090*(.182-.282)/(.090-.282)/0.182=0.2575549
296:卵の名無しさん
17/12/02 11:15:07.32 ZaK9sW49.net
# 問.
# 患者が煙草を忘れて行ったとする。
# 忘れて行った人物が女性である確率を以下のデータから計算せよ。
#
# 喫煙率
# 男性 28.2％
# 女性 9.0％
# 男女計 18.2％
P(s|m) = 0.282
P(s|f) = 0.090
P(s) = P(s|m)P(m)+P(s|f)P(f)
= P(s|m)(1-P(f)) + P(s|f)P(f)
= 0.182
P(f) = (P(s) - P(s|m))/(P(s|f) - P(s|m))
= (0.182 - 0.282)/(0.090 - 0.282)
= 0.5208333
P(f|s) = P(s|f)P(f)/P(s)
= 0.090*0.5208333/0.182
= 0.2575549
#
LR = P(s|f)|P(s|m)=0.090/0.282=0.3191489
prior.odds(f)=P(f)/(1-P(f))=0.5208333/(1-0.5208333)=1.086956
post.odds(f|s)= prior.odds(f)*LH=1.086956*0.3191489=0.3469008
P(f|s)=post.odds(f|s)/(1+post.odds(f|s))=0.3469008/(1+0.3469008)
= 0.2575548
297:卵の名無しさん
17/12/03 06:11:30.00 B6LMarvh.net
１次方程式もできないド底辺特殊シリツ医大卒の記録
URLﾘﾝｸ(imagizer.imageshack.com)
何度読んでも馬鹿すぎる。
男女別の割合と全体での割合から男女比が計算できるとも思わないとは。
なんでこんなのが大学に入れるわけよ？
裏口入学以外に説明がつく？
中学生でも解ける一次方程式の問題だろ。
それすらできない馬鹿が自信を持って発言。
>患者の男女比が必要なのもわからないのか？
だとさ。
URLﾘﾝｸ(imagizer.imageshack.com)
0.2575549
と答を書いてやったら
>単位も書かずに答えだとか…
ド底辺シリツ医大では確率に単位があるらしいぞｗｗｗ
何でこんな馬鹿が大学に入れるわけ？
裏口入学以外に説明がつく？
URLﾘﾝｸ(imagizer.imageshack.com)
298:卵の名無しさん
17/12/03 08:22:35.16 Egt6Q5KK.net
１次方程式もできないド底辺特殊シリツ医大卒の記録
URLﾘﾝｸ(imagizer.imageshack.com)
何度読んでも馬鹿すぎる。
男女別の割合と全体での割合から男女比が計算できるとも思わないとは。
なんでこんなのが大学に入れるわけよ？
裏口入学以外に説明がつく？
中学生でも解ける一次方程式の問題だろ。
シリツ医大には二次方程式が解けないやつがいると言ってた えなりかずき もビックリだろね。
それすらできない馬鹿が自信を持って発言。
>患者の男女比が必要なのもわからないのか？
だとさ。
URLﾘﾝｸ(imagizer.imageshack.com)
求める確率を
0.2575549
と答を書いてやったら
>単位も書かずに答えだとか…
ド底辺シリツ医大では確率に単位があるらしいぞｗｗｗ
何でこんな馬鹿が大学に入れるわけ？
裏口入学以外に説明がつく？
URLﾘﾝｸ(imagizer.imageshack.com)
299:卵の名無しさん
17/12/03 10:56:29.84 qW8l0b6t.net
# ある仮想の難治疾患患者25人従来薬を投与して3人治癒した。
# 新薬が登場して3人に投与したところ治癒した人はいなかった。
# この新薬を継続して使う価値があるかどうか検討せよ。
別バージョン
# 巨乳女子大で25人に声をかけたら3人が誘いにのった。
# 桃尻女子大で3人に声をかけたら誰も誘いにのらなかった。
# どちらが口説きやすいか検討せよ。
JAGSでMCMCして治癒率の確率密度関数を描くとこうなる。
URLﾘﾝｸ(i.imgur.com)
治癒率差の不偏推定量は
> mean(dif)
[1] -0.05136971
54％が負
> c(mean(dif<0),mean(dif>0))
[1] 0.5395 0.4605
５％幅の違いは同等扱いにすると
> c(mean(dif<ROPE[1]),mean(ROPE[1]<dif & dif<ROPE[2]), mean(dif>ROPE[2]))
[1] 0.4834 0.1236 0.3930
と計算できる。
98％HDIは
> HDInterval::hdi(dif)
lower upper
-0.4247349 0.2535357
と０を挟む。
RのパッケージBESTを改造して、治癒率の差の関数密度をかくと
URLﾘﾝｸ(i.imgur.com)
ゆえに新薬は無効とはいえないだけでなく、不偏推定量から従来薬を凌駕する可能性が54％ある。
300:卵の名無しさん
17/12/04 03:24:01.71 mVXTI5F+.net
# ある仮想の難治疾患患者25人従来薬を投与して3人治癒した。
# 新薬が登場して3人に投与したところ治癒した人はいなかった。
# この新薬を継続して使う価値があるかどうか検討せよ。
別バージョン
# 巨乳女子大で25人に声をかけたら3人が誘いにのった。
# 桃尻女子大で3人に声をかけたら誰も誘いにのらなかった。
# どちらが口説きやすいか検討せよ。
JAGSでMCMCして治癒率の確率密度関数を描くとこうなる。
URLﾘﾝｸ(i.imgur.com)
治癒率差の不偏推定量は
> mean(dif)
[1] -0.05136971
54％が負
> c(mean(dif<0),mean(dif>0))
[1] 0.5395 0.4605
５％幅の違いは同等扱いにすると
> c(mean(dif<ROPE[1]),mean(ROPE[1]<dif & dif<ROPE[2]), mean(dif>ROPE[2]))
[1] 0.4834 0.1236 0.3930
と計算できる。
98％HDIは
> HDInterval::hdi(dif)
lower upper
-0.4247349 0.2535357
と０を挟む。
RのパッケージBESTを改造して、治癒率の差の確率密度をかくと
URLﾘﾝｸ(i.imgur.com)
ゆえに新薬は無効とはいえないだけでなく、不偏推定量から従来薬を凌駕する可能性が54％ある。
301:卵の名無しさん
17/12/04 12:27:18.87 dllejky7.net
# ある仮想の難治疾患患者25人に従来薬を投与して3人治癒した。
# 新薬が登場して3人に投与したところ治癒した人はいなかった。
# この新薬を継続して使う価値があるかどうか検討せよ。
別バージョン
# 巨乳女子大で25人に声をかけたら3人が誘いにのった。
# 桃尻女子大で3人に声をかけたら誰も誘いにのらなかった。
# どちらが口説きやすいか検討せよ。
JAGSでMCMCして治癒率の確率密度関数を描くとこうなる。
URLﾘﾝｸ(i.imgur.com)
治癒率差の不偏推定量は
> mean(dif)
[1] -0.05136971
54％が負
> c(mean(dif<0),mean(dif>0))
[1] 0.5395 0.4605
５％幅の違いは同等扱いにすると
> c(mean(dif<ROPE[1]),mean(ROPE[1]<dif & dif<ROPE[2]), mean(dif>ROPE[2]))
[1] 0.4834 0.1236 0.3930
と計算できる。
95％HDIは
> HDInterval::hdi(dif)
lower upper
-0.4247349 0.2535357
と０を挟む。
RのパッケージBESTを改造して、治癒率の差の確率密度分布をかくと
URLﾘﾝｸ(i.imgur.com)
ゆえに新薬は無効とはいえないだけでなく、不偏推定量から従来薬を凌駕する可能性が54％ある。
302:卵の名無しさん
17/12/08 20:35:21.62 raU+TCc7.net
In summary, when there is interaction, then the influence of the individual predictors can not be summarized by their individual regression coefficients alone, because those coefficients only describe the influence when the other variables are at zero.
A careful analyst considers credible slopes across a variety of values for the other predictors.
Notice that this is true even though the interaction coefficient did not exclude zero from its 95% HDI.
In other words, if you include an interaction term, you cannot ignore it even if its marginal posterior distribution includes zero.
303:卵の名無しさん
17/12/10 13:21:50.45 xzT2/Bky.net
seqn<-function(n=5,N=100,p=0.5){ # N回のうちn回以上続けて表がでるか？
rn=rbinom(N,1,p)
count=0
for(i in 1:N){
if(rn[i] & count<n){
count=count+1
}
else{
if(count==n) {return(TRUE)}
else{
count=0
}
}
}
return(count==n)
}
mean(replicate(10^5,seqn()))
f <- function(n) mean(replicate(10^4,seqn(n)))
nn=2:20
yy=sapply(nn,f)
plot(nn,yy,pch=19,xlab='sequential heads',ylab='Proportion')
abline(h=0.05,lty=3)
f(9)
f(10)
304:卵の名無しさん
17/12/10 14:45:43.48 xzT2/Bky.net
# 最頻値M 平均m 分散v のガンマ分布を作る
Mv2sr <- function(M,v){
shape=(M^2 +2*v+sqrt(M^2*(M^2+4*v)))/(2*v)
rate= (M^2+ sqrt(M^2*(M^2+4*v)))/(2*M*v)
c(shape=shape,rate=rate)
}
Mv2sr(1,1)
sr2mMv <- function(shape,rate){
c(mean=shape/rate,mode=(shape-1)/rate,var=shape/(rate^2))
}
sr2mMv(2.618,1.618)
mv2sr <- function(mean,var){
rate=mean/var
shape=mean*rate
c(shape=shape,rate=rate)
}
mv2sr(1.618,1)
305:卵の名無しさん
17/12/12 06:01:50.72 6Uyksjmd.net
URLﾘﾝｸ(imagizer.imageshack.com)
306:卵の名無しさん
17/12/12 20:18:02.76 6Uyksjmd.net
f <- function(i){
re=i+0:(k-1)
re=re%%n
re[which(re==0)]=n
return(re)
}
g <- function(x) (x+1)%%2
h <- function(i,b){
idx=f(i)
b[idx]=g(b[idx])
return(b)
}
i <- function(v){
tmp=a
for(w in v){
tmp=h(w,tmp)
}
return(tmp)
}
307:卵の名無しさん
17/12/12 20:19:54.91 6Uyksjmd.net
n=7
k=3
a=rep(0,7) #7枚全部裏のとき
f <- function(i){
re=i+0:(k-1)
re=re%%n
re[which(re==0)]=n
return(re)
}
g <- function(x) (x+1)%%2
h <- function(i,b){
idx=f(i)
b[idx]=g(b[idx])
return(b)
}
i <- function(v){
tmp=a
for(w in v){
tmp=h(w,tmp)
}
return(tmp)
}
308:卵の名無しさん
17/12/12 20:21:21.43 6Uyksjmd.net
sc6=t(combn(n,6))
sc6p=numeric(nrow(sc6))
for(j in 1:nrow(sc6)){
sc6p[j]=prod(i(sc6[j,]))
}
any(sc6p==1) #6回でも無理
sc7=t(combn(n,7))
sc7p=numeric(nrow(sc7))
for(j in 1:nrow(sc7)){
sc7p[j]=prod(i(sc7[j,]))
}
any(sc7p==1) # TRUE! 7回で全部表にできる
sc7[which(sc7p==1),]
実行してみる(0:裏 1:表)
0000000
1110000
1001000
1010100
1011010
1011101
0011110
1111111
309:卵の名無しさん
17/12/16 15:56:30.23 LibxCzfo.net
平均・標準偏差が以下の４群の多重比較
> tapply(Y,Group,mean)
A B C D
97 99 102 104
> tapply(Y,Group,sd)
A B C D
8 1 1 8
URLﾘﾝｸ(i.imgur.com)
URLﾘﾝｸ(i.imgur.com)
> kruskal.test(Y~Group)
Kruskal-Wallis rank sum test
data: Y by Group
Kruskal-Wallis chi-squared = 25.325, df = 3, p-value = 0.0000132
> pairwise.t.test(Y,Group,p.adjust.method = 'holm', pool.sd = FALSE)
Pairwise comparisons using t tests with non-pooled SD
data: Y and Group
A B C
B 0.472 - -
C 0.023 0.00000000000071 -
D 0.020 0.023 0.472
P value adjustment method: holm
310:卵の名無しさん
17/12/16 15:58:37.94 LibxCzfo.net
> pairwise.t.test(Y,Group,p.adjust.method = 'bon', pool.sd = FALSE)
Pairwise comparisons using t tests with non-pooled SD
data: Y and Group
A B C
B 1.000 - -
C 0.034 0.00000000000071 -
D 0.024 0.034 1.000
P value adjustment method: bonferroni
> pairwise.t.test(Y,Group,p.adjust.method = 'fdr', pool.sd = FALSE)
Pairwise comparisons using t tests with non-pooled SD
data: Y and Group
A B C
B 0.2362 - -
C 0.0086 0.00000000000071 -
D 0.0086 0.0086 0.2362
P value adjustment method: fdr
どの補正でも有意差ありは、A-C,A-D,B-C,B-D
有意差なしは A-B,C-D
311:卵の名無しさん
17/12/16 16:03:58.73 LibxCzfo.net
等分散でないのでWilcoxは参考程度だが、結果は同じ。
> pairwise.wilcox.test(Y,Group,p.ad='holm')
Pairwise comparisons using Wilcoxon rank sum test
data: Y and Group
A B C
B 0.231 - -
C 0.005 0.000000000073 -
D 0.019 0.062 0.533
P value adjustment method:
312:holm
313:卵の名無しさん
17/12/16 16:06:39.96 LibxCzfo.net
ｐ値で考えないから補正も不要のベイズの方が直感的だな。
URLﾘﾝｸ(i.imgur.com)
314:卵の名無しさん
17/12/16 16:24:11.25 LibxCzfo.net
分散が大きく異なり、有意差がないときはmcmcでの図示で理解が深まる。
URLﾘﾝｸ(i.imgur.com)
315:卵の名無しさん
17/12/16 16:51:04.28 LibxCzfo.net
プールした標準偏差をつかうとB-Cが有意でなくなる。
> pairwise.t.test(Y,Group,p.adjust.method = 'holm',pool.sd = TRUE)
Pairwise comparisons using t tests with pooled SD
data: Y and Group
A B C
B 0.4547 - -
C 0.0155 0.2147 -
D 0.0003 0.0155 0.4547
P value adjustment method: holm
MCMCのモデルで標準偏差をグループ共通にすると95％HDIが０を跨いでしまう。
URLﾘﾝｸ(i.imgur.com)
316:卵の名無しさん
17/12/16 16:52:42.18 LibxCzfo.net
グループごとに固有の標準偏差を持つとすると
URLﾘﾝｸ(i.imgur.com)
で HDIが０を跨がない。
317:卵の名無しさん
17/12/16 17:05:01.28 LibxCzfo.net
pooled.sdモデル
URLﾘﾝｸ(i.imgur.com)
固有sdモデル
URLﾘﾝｸ(i.imgur.com)
どちらが現実をよりよく説明するかということだな。
318:卵の名無しさん
17/12/16 23:23:05.01 LibxCzfo.net
# 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)
}
319:卵の名無しさん
17/12/17 04:09:30.35 o6qOChW2.net
# N個のクジでr個めで初めてあたった時のN個内の当たり数の推測
N=100 ; r=5
pmf <- function(x) (1-x/N)^(r-1)*x/N # dnbinom(r-1,1,x/N) ; dgeom(r-1,x/N)
curve((1-x/N)^(r-1)*x/N,0,N)
AUC=integrate(pmf,0,N)$value
pdf <- function(x) pmf(x)/AUC
source('tools.R')
pdf2hdi(pdf,0,N)
320:卵の名無しさん
17/12/17 13:39:11.48 o6qOChW2.net
>>282
female_propotion=100/(192)
kappa=10
(AB=betaABfromMeanKappa(female_propotion,kappa))
alpha=AB$a ; beta=AB$b
# alpha=1 ; beta=1
curve(dbeta(x,alpha,beta),xlab='female/male ratio')
data=list(
smoker_female=0.090,
smoker_male=0.282,
smoker=0.182,
alpha=alpha,
beta=beta
)
321:卵の名無しさん
17/12/17 13:39:44.90 o6qOChW2.net
stanStrings='
data{
real<lower=0,upper=1> smoker_female;
real<lower=0,upper=1> smoker_male;
real<lower=0,upper=1> smoker;
real<lower=0> alpha;
real<lower=0> beta;
}
parameters{
real<lower=0,upper=1> female;
}
model{
female ~ beta(alpha,beta);
}
generated quantities{
real<lower=0,upper=1> female_smoker;
female_smoker = smoker_female*female/smoker;
}
'
stanModel=stan_model(model_code = stanStrings)
fit=sampling(stanModel,data=data,seed=123,iter=50000,warmup=5000,chains=4)
print(fit,digits=4)
322:卵の名無しさん
17/12/17 18:48:25.00 o6qOChW2.net
テキストで解説したあるグラフが自分で再現できないと気になるね。
ようやく完成。
URLﾘﾝｸ(i.imgur.com)
べつに分布を９０度回転させて表示させなくてもいいのだが。
323:卵の名無しさん
17/12/17 20:32:29.49 o6qOChW2.net
http//i.imgur.com/fzzGWoz.png
324:卵の名無しさん
17/12/18 08:45:56.18 51j+AsC2.net
URLﾘﾝｸ(i.imgur.com)
325:卵の名無しさん
17/12/18 22:32:27.86 v8MzYeil.net
# 薬剤ｙは3人目で治癒、薬剤ｇは10人中2人が治癒、どちらの有効性が高いか？
stanStrings='
data{
int r; //3
int z; //3
int N; //10
int<lower=0,upper=1> y[r]; //c(0,0,1)
}
parameters{
real<lower=0,upper=1>p[2]; //p[1]:yuruyuru, p[2]:gabagaba
}
model{
y ~ bernoulli(p[1]);
z ~ binomial(N,p[2]);
}
generated quantities{
real<lower=0,upper=1> yuru;
real<lower=0,upper=1> gaba;
real diff;
yuru = (1-p[1])^(r-1)*p[1];
gaba = choose(N,z)*p[2]^z*(1-p[2])^(N-z)
326:; diff = p[1]-p[2]; } ' data=list(r=3,z=3,N=10,y=c(0,0,1)) stanmodel=stan_model(model_code = stanStrings) fit=sampling(stanmodel,data=data,seed=123) print(fit,digits=4)
327:卵の名無しさん
17/12/18 23:42:11.61 v8MzYeil.net
>>310
声を掛けたら、ゆるゆる女子大r人めで開脚、がばがば女子大N人中ｚ人開脚、どっちが開脚が容易か？
という問題にした方が興味をひくなｗ
stanでの結果は これ
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
p[1] 0.399 0.001 0.200 0.066 0.241 0.386 0.542 0.802 31981 1
p[2] 0.334 0.001 0.130 0.110 0.238 0.324 0.421 0.609 32098 1
yuru 0.114 0.000 0.035 0.028 0.094 0.127 0.143 0.148 19842 1
gaba 0.195 0.000 0.071 0.029 0.150 0.218 0.255 0.267 20013 1
diff 0.065 0.001 0.239 -0.377 -0.106 0.058 0.231 0.541 32024 1
lp__ -12.079 0.008 1.063 -14.941 -12.503 -11.751 -11.313 -11.031 16102 1
p[1]:ゆるゆる女子大開脚率
p[2]:がばがば女子大開脚率
diff = p[1]-p[2];
328:卵の名無しさん
17/12/19 13:22:41.16 t02o1U8G.net
>>311
# ゆるゆる女子大生 r 人めではじめて開脚、がばがば女子大生 N 人中 ｚ 人が開脚、どっちが開脚が容易か？
r=3
z=3
N=9
とサンプルでの比率が同じとき母集団の推定平均値に差があるだろうか？
stanの出力をグラフにしてみた。
平均値で４％ほどの差が推定された。
URLﾘﾝｸ(i.imgur.com)
329:卵の名無しさん
17/12/19 17:05:34.32 t02o1U8G.net
# dnbinom(10,5,0.5) 5回表をだすまでに10回裏がでる確率
dnbinom(24-7,7,7/24)
N24=24
z7=7
nn=0:50
pp=dnbinom(nn,z7,z7/(nn+z7))
plot(z7/(nn+z7),pp,type='h',col='blue')
points(z7/N24,0, pch='+', cex=2,col=2)
330:卵の名無しさん
17/12/20 19:30:19.77 hlUsptvw.net
あるド底辺シリツ医大で入学者の裏口入学者と学力考査合格入学者の比率は１であるという帰無仮説を検定する課題が
花子と太郎に課された。
花子は５０人を調査できたら終了としてド底辺シリツ医大入学者を５０人をみつけて１８人が裏口であるという結果を得た。
帰無仮説のもとで
５０人中１８人が裏口である確率は 0.01603475
これ以下になるのは５０人中０?１８人と３２?５０人が裏口の場合なので
両側検定して
> sum(dbinom(c(0:18,32:50),50,0.5))
[1] 0.06490865
> binom.test(18,50,0.5)$p.value
[1] 0.06490865
で帰無仮説は棄却できないと結論した。
URLﾘﾝｸ(i.imgur.com)
一方、本番と十八番が好きな太郎は一人ずつ調べて１８人めの裏口がみつかったところで調査を終えることにした。
１８人めがみつかったのは花子と同じく５０人めであった。
帰無仮説のもとで
１８人がみつかるのが５０人めである確率は0.005772512
これ以下になるのは２３人以下５０人以上番目で裏口１８人めがみつかった場合なので
両側検定して
pnb=dnbinom(0:999,18,0.5)
> 1 - sum(pnb[-which(pnb<=dnbinom(50-18,18,0.5))]) # < 0.05
[1] 0.02750309
URLﾘﾝｸ(i.imgur.com)
で帰無仮説は棄却される。
331:卵の名無しさん
17/12/20 21:16:42.01 e+3oE/TR.net
コインが続けて５回裏がでたときにこのコインはイカサマコインといえるか？
５％のばらつきはイカサマとみなさないとする。
ROPE=c(0.5,0.5*1.05)
curve(dbeta(x,1+zi,1+Ni-zi))
abline(v=ROPE[1],col='gray',lty=3) ;abline(v=ROPE[2],col='gray',lty=3)
pbeta(ROPE[2],1+zi,1+Ni-zi)-pbeta(ROPE[1],1+zi,1+Ni-zi)
HDInterval::hdi(qbeta,shape1=1+zi,shape2=1+Ni-zi)
332:卵の名無しさん
17/12/20 23:11:31.75 e+3oE/TR.net
> require('HDI')
> f=function(n) HDInterval::hdi(qbeta,shape1=n+1,shape2=1)[1]
> f(5)
lower
0.6069622
> f(10)
lower
0.7615958
> re=sapply (1:100,f)
> names (re)=NULL
> re
[1] 0.2236068 0.3684031 0.4728708 0.5492803 0.6069622 0.6518363 0.6876560
[8] 0.7168712 0.7411344 0.7615958 0.7790778 0.7941833 0.8073638 0.8189637
[15] 0.8292503 0.8384339 0.8466824 0.8541315 0.8608917 0.8670541 0.8726946
[22] 0.8778766 0.8826538 0.8870719 0.8911696 0.8949808 0.8985343 0.9018554
[29] 0.9049661 0.9078859 0.9106318 0.9132188 0.9156604 0.9179684 0.9201535
[36] 0.9222253 0.9241923 0.9260624 0.9278425 0.9295389 0.9311574 0.9327032
[43] 0.9341812 0.9355957 0.9369507 0.9382499 0.9394966 0.9406940 0.9418449
[50] 0.9429520 0.9440178 0.9450445 0.9460342 0.9469889 0.9479105 0.9488005
[57] 0.9496607 0.9504924 0.9512971 0.9520760 0.9528305 0.9535615 0.9542703
[64] 0.9549577 0.9556248 0.9562724 0.9569014 0.9575126 0.9581067 0.9586843
[71] 0.9592463 0.9597932 0.9603256 0.9608441 0.9613492 0.9618415 0.9623214
[78] 0.9627893 0.9632458 0.9636912 0.9641260 0.9645504 0.9649650 0.9653699
[85] 0.9657656 0.9661524 0.9665305 0.9669003 0.9672620 0.9676158 0.9679621
[92] 0.9683011 0.9686330 0.9689580 0.9692763 0.9695882 0.9698938 0.9701933
[99] 0.9704869 0.9707748
>
333:卵の名無しさん
17/12/20 23:30:49.76 e+3oE/TR.net
> g=function (n) qbeta(0.05,shape1=n+1,shape2=1)
> g(5)
[1] 0.6069622
> g(10)
[1] 0.7615958
> sapply(1:100,g)
[1] 0.2236068 0.3684031 0.4728708 0.5492803 0.6069622 0.6518363 0.6876560
[8] 0.7168712 0.7411344 0.7615958 0.7790778 0.7941833 0.8073638 0.8189637
[15] 0.8292503 0.8384339 0.8466824 0.8541315 0.8608917 0.8670541 0.8726946
[22] 0.8778766 0.8826538 0.8870719 0.8911696 0.8949808 0.8985343 0.9018554
[29] 0.9049661 0.9078859 0.9106318 0.9132188 0.9156604 0.9179684 0.9201535
[36] 0.9222253 0.9241923 0.9260624 0.9278425 0.9295389 0.9311574 0.9327032
[43] 0.9341812 0.9355957 0.9369507 0.9382499 0.9394966 0.9406940 0.9418449
[50] 0.9429520 0.9440178 0.9450445 0.9460342 0.9469889 0.9479105 0.9488005
[57] 0.9496607 0.9504924 0.9512971 0.9520760 0.9528305 0.9535615 0.9542703
[64] 0.9549577 0.9556248 0.9562724 0.9569014 0.9575126 0.9581067 0.9586843
[71] 0.9592463 0.9597932 0.9603256 0.9608441 0.9613492 0.9618415 0.9623214
[78] 0.9627893 0.9632458 0.9636912 0.9641260 0.9645504 0.9649650 0.9653699
[85] 0.9657656 0.9661524 0.9665305 0.9669003 0.9672620 0.9676158 0.9679621
[92] 0.9683011 0.9686330 0.9689580 0.9692763 0.9695882 0.9698938 0.9701933
[99] 0.9704870 0.9707748
>
334:卵の名無しさん
17/12/21 02:15:57.98 NCbCbV7K.net
> binom::binom.confint(18,50)
method x n mean lower upper
1 agresti-coull 18 50 0.3600000 0.2410278 0.4989496
2 asymptotic 18 50 0.3600000 0.2269532 0.4930468
3 bayes 18 50 0.3627451 0.2343802 0.4940800
4 cloglog 18 50 0.3600000 0.2306356 0.4908871
5 exact 18 50 0.3600000 0.2291571 0.5080686
6 logit 18 50 0.3600000 0.2399736 0.5005239
7 probit 18 50 0.3600000 0.2375867 0.4988707
8 profile 18 50 0.3600000 0.2363864 0.4976324
9 lrt 18 50 0.3600000 0.2363786 0.4976328
10 prop.test 18 50 0.3600000 0.2328502 0.5085700
11 wilson 18 50 0.3600000 0.2413875 0.4985898
> HDInterval::hdi(qbeta,shape1=19,shape2=33)
lower upper
0.2379677 0.4956588
attr(,"credMass")
[1] 0.95
>
335:卵の名無しさん
17/12/21 05:06:47.13 NCbCbV7K.net
> 0.05^(1/(1:100))
[1] 0.0500000 0.2236068 0.3684031 0.4728708 0.5492803 0.6069622 0.6518363
[8] 0.6876560 0.7168712 0.7411344 0.7615958 0.7790778 0.7941833 0.8073638
[15] 0.8189637 0.8292503 0.8384339 0.8466824 0.8541315 0.8608917 0.8670541
[22] 0.8726946 0.8778766 0.8826538 0.8870719 0.8911696 0.8949808 0.8985343
[29] 0.9018554 0.9049661 0.9078859 0.9106318 0.9132188 0.9156604 0.9179684
[36] 0.9201535 0.9222253 0.9241923 0.9260624 0.9278425 0.9295389 0.9311574
[43] 0.9327032 0.9341812 0.9355957 0.9369507 0.9382499 0.9394966 0.9406940
[50] 0.9418449 0.9429520 0.9440178 0.9450445 0.9460342 0.9469889 0.9479105
[57] 0.9488005 0.9496607 0.9504924 0.9512971 0.9520760 0.9528305 0.9535615
[64] 0.9542703 0.9549577 0.9556248 0.9562724 0.9569014 0.9575126 0.9581067
[71] 0.9586843 0.9592463 0.9597932 0.9603256 0.9608441 0.9613492 0.9618415
[78] 0.9623214 0.9627893 0.9632458 0.9636912 0.9641260 0.9645504 0.9649650
[85] 0.9653699 0.9657656 0.9661524 0.9665305 0.9669003 0.9672620 0.9676158
[92] 0.9679621 0.9683011 0.9686330 0.9689580 0.9692763 0.9695882 0.9698938
[99] 0.9701933 0.9704870
336:卵の名無しさん
17/12/21 05:20:38.89 NCbCbV7K.net
> 0.05^(1/(1:100+1))
[1] 0.2236068 0.3684031 0.4728708 0.5492803 0.6069622 0.6518363 0.6876560
[8] 0.7168712 0.7411344 0.7615958 0.7790778 0.7941833 0.8073638 0.8189637
[15] 0.8292503 0.8384339 0.8466824 0.8541315 0.8608917 0.8670541 0.8726946
[22] 0.8778766 0.8826538 0.8870719 0.8911696 0.8949808 0.8985343 0.9018554
[29] 0.9049661 0.9078859 0.9106318 0.9132188 0.9156604 0.9179684 0.9201535
[36] 0.9222253 0.9241923 0.9260624 0.9278425 0.9295389 0.9311574 0.9327032
[43] 0.9341812 0.9355957 0.9369507 0.9382499 0.9394966 0.9406940 0.9418449
[50] 0.9429520 0.9440178 0.9450445 0.9460342 0.9469889 0.9479105 0.9488005
[57] 0.9496607 0.9504924 0.9512971 0.9520760 0.9528305 0.9535615 0.9542703
[64] 0.9549577 0.9556248 0.9562724 0.9569014 0.9575126 0.9581067 0.9586843
[71] 0.9592463 0.9597932 0.9603256 0.9608441 0.9613492 0.9618415 0.9623214
[78] 0.9627893 0.9632458 0.9636912 0.9641260 0.9645504 0.9649650 0.9653699
[85] 0.9657656 0.9661524 0.9665305 0.9669003 0.9672620 0.9676158 0.9679621
[92] 0.9683011 0.9686330 0.9689580 0.9692763 0.9695882 0.9698938 0.9701933
[99] 0.9704870 0.9707748
>
337:卵の名無しさん
17/12/22 02:11:19.86 J9UAx7pH.net
シンプソンのパラドックス
ある仮想疾患の治癒率
軽症 重症
国立大学 10/10 10/90
底辺私立 70/90 0/10
自然経過 40/50 5/50
国立大学の方が軽症・重症とも成績がよいが
総数比較では底辺私立の方が成績がよい。
この疾患は自然治癒率が45%とされています。
この疾患の底辺私立での治癒率は70%です。
これに�
338:ﾎして国立大学での治癒率はわずか20%です。 という記述も嘘ではないね
339:卵の名無しさん
17/12/22 12:38:20.39 FPEZRkaT.net
f <- function(n=10,alpha=1,beta=1,Print=FALSE){
N=n
z=n
if(Print) {
bayes=binom::binom.bayes(z,N, prior.shape1=alpha,prior.shape2=beta)
show(binom::binom.bayes.densityplot(bayes))
}
hdi=HDInterval::hdi(qbeta, shape1=alpha+z, shape2=beta+N-z)
return (c(lower=hdi[1],mean=(alpha+z)/(alpha+beta+N),
mode=(alpha+z-1)/(alpha+beta+N-2), upper=hdi[2]))
}
f(10,P=TRUE)
nn=1:30
yy=sapply(nn,function(x)f(x,Print=FALSE)[1])
plot(nn,yy,pch=19,xlab='裏口バカ連続合格数',ylab='裏口確率信頼区間下限')
curve((0.05)^(1/(x+1)),add=TRUE,lty=3) # 0.05の（合格者数＋１）乗根
340:卵の名無しさん
17/12/22 20:41:25.07 JBU22EfC.net
# N回続けて裏、事前分布はmode値0.5, 集中度(形状母数和)=kappa
source('tools.R')
N=5 ; z=5
Bayes(N,z,alpha=1,beta=1,Print=TRUE) ; 0.05^(1/(N+1)) # N=z
# 事前分布が最頻値0.5で集中度(κ=α+β)のとき事後分布の関係
.mode=0.5
Kappa2Bayes <- function(kappa,.mode=0.5){
AB=betaABfromModeKappa(.mode,kappa)
Bayes(N,z,alpha=AB[[1]],beta=AB[[2]])
}
K=seq(2,1000,by=0.5)
K=K[-1]
res=sapply(K,Kappa2Bayes)
Mat=as.matrix(res)
plot(K,Mat['lower',],type='l',xlab=bquote(kappa),
ylab='Probability', ylim=c(0,1),lty=3)
lines(K,Mat['mean',],col=4,lwd=4)
lines(K,Mat['mode',],col=2,lwd=2)
lines(K,Mat['upper',],lty=3)
legend('bottomright',bty='n',legend=c('mean','mode','upr','lwr'),
col=c(4,2,1,1),lty=c(1,1,3,3),lwd=c(4,2,1,1))
341:卵の名無しさん
17/12/22 21:31:09.94 JBU22EfC.net
N(=100)回コインをなげてｎ(=5回)以上続けて表がでる
seqn<-function(n=5,N=100,p=0.5){
rn=rbinom(N,1,p)
count=0
for(i in 1:N){
if(rn[i] & count<n){
count=count+1
}
else{
if(count==n) {return(TRUE)}
else{
count=0
}
}
}
return(count==n)
}
mean(replicate(10^6,seqn()))
> mean(replicate(10^6,seqn()))
[1] 0.810223
342:卵の名無しさん
17/12/23 05:59:24.55 RXgTWS3Z.net
pooledVariance <- function(...) {
args = list(...)
n.args=length(args)
ss2=0
df=0
for(i in 1:n.args){
ss2 = ss2 + var(args[[i]])*(length(args[[i]])-1)
df = df + (length(args[[i]])-1)
}
ss2/df
}
effectsize <- function(y1,y2){
diff=mean(y1)-mean(y2)
var=(var(x1)*(length(x1)-1)+ var(x2)*(length(x2)-1))/(length(c(y1,y2))-2)
sd=sqrt(var)
diff/sd
}
library(effsize)
cohen.d()
343:卵の名無しさん
17/12/23 06:36:30.75 RXgTWS3Z.net
solve[{M=(a-1)/(a+b-2), V=a*b/((a+b)^2*(a+b+1))},{a,b}]
344:卵の名無しさん
17/12/23 17:39:59.62 K1KZza+Z.net
> 1/(1-(1-0.99)^(1/317))
[1] 69.33689
1 / ( 1- n√(1-confidence.level) )
345:卵の名無しさん
17/12/23 18:01:15.54 K1KZza+Z.net
# 1 / ( 1- n√(1-confidence.level) )
confidence.level=0.95
rule3 <- function(n,confidence.level=0.95){ # n人に1人の副作用
p=1/n
q=1-p # q^n.sample < 1-confidence.level
n.sample = log(1-confidence.level)/log(q)
return(n.sample)
}
nn=seq(1,10000,by=100)
plot(nn,sapply(nn,rule3))
curve(3*x,col=2,add=TRUE)
plot(nn,sapply(nn,function(x) rule3(x,conf=0.99)))
lm( sapply(nn,function(x) rule3(x,conf=0.99)) ~ nn + 0)
curve(4.605*x,col=4,add=TRUE)
346:卵の名無しさん
17/12/24 20:11:33.25 SYFul+nD.net
サイコロの１の目が何回続けてでたらイカサマサイコロかを考えて暇つぶし。
このスレの趣旨ｗに合わせてこんな問題にしてみた。
ド底辺学力学生がド底辺特殊シリツ医大を受験したとする。
試験は五者択一で三教科150問、合格ラインは６割とされる。
中学数学すらできないド底辺学力ゆえ正解できるのは偶然に頼るしか�
347:ﾈく、 正答率は概ね1/５で、その日の運で1/4から1/6と推定されている。 これを、正答確率は最頻値1/５で1/6から1/4の間に正答する確率の95%があると設定する。 このド底辺が合格したとする。150問中90問以上正答したことになる。 これがイカサマ入試である確率を求めよ。 事前確率のβ分布のパラメータ算出がやや手間だが、あとはルーティン作業 JAGSを使ってMCMCでの結果 http://i.imgur.com/V7TaBG7.png 解析解： 0.9994608 とほぼ一致。
348:卵の名無しさん
17/12/25 08:52:06.92 Nj//P1mP.net
require(rjags)
N=10
z=0
y=c(rep(1,z),rep(0,N-z))
ph=c(1,2/3,1/2,1/5)
pc=c(2/100,50/100,40/100,8/100)
names(ph)=c('特待生','学力合格','加点合格','ガチの裏')
names(pc)=c('特待生','学力合格','加点合格','ガチの裏')
dataList5=list(N=N,y=y,ph=ph,pc=pc)
# JAGS model
modelString5 ="
model {
for(i in 1:N){
y[i] ~ dbern(ph[coin])
}
coin ~ dcat(pc[])
}
"
writeLines(modelString5,'TEMPmodel.txt')
jagsModel5=jags.model('TEMPmodel.txt',data=dataList5)
codaSamples5=coda.samples(jagsModel5,var=c('coin'),n.iter=100000,na.rm=TRUE)
summary(codaSamples5)
js5=as.matrix(codaSamples5)
re=numeric()
for(i in 1:4) re[i]=mean(js5==i)
dat=data.frame(割合=round(re*100,3))
rownames(dat)=names(pc)
dat
349:卵の名無しさん
17/12/25 16:31:55.67 UMwuImpO.net
頻度主義統計の謎。
立方体からなるサイコロの目のでる確率はすべて等しく1/6である、を帰無仮説とする。
そのサイコロをふって1の目がでた。2回目は2の目がでた。
その確率は1/6*1/6で1/36=0.02778 < 0.05だから帰無仮説は棄却される。
どの目の組合せでも同じく帰無仮説は棄却される。
頻度主義統計のもとではすべてのサイコロはいびつである。
350:卵の名無しさん
17/12/27 10:24:07.92 un/eaZi1.net
a005=0.05
n100=100
HDInterval::hdi(qbeta,shape1=n100+.a,shape2=.b)
qbeta(a005,n100+.a,.b)
a005=0.05
# p^n100 < a005
# p < a005^(1/n100)
GolgoLowerLimit <- function(a005,n100=100){ # Golgo lower limit
c(a005^(1/(n100+1)),qbeta(a005,n100+1,1))
}
GolgoLowerLimit(0.05)
AHO <- function(a00,n100=100,shoot=10000){
(a005^(1/n100)+1)/2*shoot
}
AHO(0.05)
x=seq(0.001,0.1,by=0.001)
plot(x,sapply(x,function(x) AHO(x,100,10000)),type='l',lwd=2,
las=1,ylab='「平均値」',xlab='危険率')
351:卵の名無しさん
17/12/27 10:57:46.05 un/eaZi1.net
data{
int Npip; // 6
real alpha[Npip]; // c(1,1,1,1,1,1)
int Ntotal; // length(y)
int y[Ntotal];
}
parameters{
simplex[Npip] pi;
}
model{
for(i in 1:Ntotal){
y[i] ~ categorical(pi);
pi ~ dirichlet(alpha)
}
}
352:卵の名無しさん
17/12/28 15:27:03.23 t4TEzXKz.net
source('tools.R')
RBI <- function(a=1,b=1,zi=1,Ni=1,zj=0,Nj=2,ROPE=NULL,Print=TRUE){
(y=c(rep(1,zi),rep(0,Ni-zi),rep(1,zj),rep(0,Nj-zj)))
(s=as.numeric(factor(c(rep('D',Ni),rep('U',Nj)))))
myData=data.frame(y=y,s=s)
Ntotal = length(y)
Nsubj = length(unique(s))
dataList = list(
y = y ,
s = s ,
Ntotal = Ntotal ,
Nsubj = Nsubj
)
# JAGS model
modelString = paste0("
model {
for ( i in 1:Ntotal ) {
y[i] ~ dbern( theta[s[i]] )
}
for ( sIdx in 1:Nsubj ) {
theta[sIdx] ~ dbeta(", a,',' , b," )
}
}
")
# close modelString
353:卵の名無しさん
17/12/28 15:41:40.10 t4TEzXKz.net
writeLines( modelString , con="TEMPmodel.txt" )
jagsModel = jags.model( file="TEMPmodel.txt" , data=dataList)
codaSamples = coda.samples( jagsModel , variable="theta", n.iter=10000 )
mcmcMat=as.matrix(
354:codaSamples) if(Print){ par(mfrow=c(2,1)) curve(dbeta(x,a,b),xlab=paste0('Beta(',a,',',b,')'),bty='n',yaxt='n',ylab='', type='h', col='blue') print(plotPost2(mcmcMat[,1]-mcmcMat[,2],compVal=0,ROPE=ROPE, cenT='mean',xlab=bquote(Delta),cex=1,showCurve=FALSE)) } dif=mcmcMat[,1]-mcmcMat[,2] print(c(HDInterval::hdi(dif),mean=mean(dif))) invisible(dif) }
355:卵の名無しさん
17/12/29 18:23:01.86 lWMSqtax.net
URLﾘﾝｸ(image.slidesharecdn.com)
356:卵の名無しさん
17/12/29 18:23:24.53 lWMSqtax.net
URLﾘﾝｸ(image.slidesharecdn.com)
357:卵の名無しさん
17/12/30 02:13:47.38 8e/jZDFc.net
BF01<- function(n,z,p0,a=1,b=1,p1=0.5) {
bf=gamma(a)*gamma(b)/gamma(a+b) *
gamma(a+b+n)/gamma(a+z)/gamma(b+n-z) *
p0^z*(1-p0)^(n-z)
pri=p1/(1-p1)
pos=pri*bf
post.prob=pos/(1+pos)
c(BayesFactor=bf, PostProb=post.prob)
}
BF01(5,4,0.5,1,1)
BF01(7,7,0.5,1,1)
358:卵の名無しさん
17/12/31 08:13:46.23 2OFZ1/Lf.net
URLﾘﾝｸ(to-kei.net)
359:卵の名無しさん
17/12/31 08:29:05.13 2OFZ1/Lf.net
a=1
b=1
n=5
z=5
p=0.5
dataList=list(
a=a,
b=b,
n=n,
z=z,
p=p
)
modelString,='
model{
bf=gamma(a)*gamma(b)/gamma(a+b) *
gamma(a+b+n)/gamma(a+z)/gamma(b+n-z) *
p^z*(1-p)^(n-z) # bf= (n+1)*p^n
pri=p0/(1-p0)
pos=pri*bf
post_prob=pos/(1+pos)
p0 ~ beta(1,1)
}
'
360:卵の名無しさん
17/12/31 09:49:21.42 2OFZ1/Lf.net
a=1
b=1
n=5
z=n
p=0.5
dataList=list(
a=a,
b=b,
n=n,
z=z,
p=p
)
modelString,='
model{
bf= (n+1)*p^n
pri=p0/(1-p0)
pos=pri*bf
post_prob=pos/(1+pos)
p0 ~ beta(1,1)
}
'
writeLines(modelString,'TEMPmodel.txt')
jagsModel=jags.model('TEMPmodel.txt', data=dataList)
codaSamples=coda.samples(jagsModel,n.iter=10000,var=c
('p0'))
js=as.matrix(codaSamples)
361:卵の名無しさん
17/12/31 10:18:36.31 2OFZ1/Lf.net
５試合連続で勝敗予想的中なら頻度主義では予知能力あるとされる。p=0.03125 < 0.05
URLﾘﾝｸ(to-kei.net)
ベイズでやってみるなら
的中率が1/2である確率は一応分布に従う（事前分布）として
５試合連続的中した後の的中率事後分布がどうなるかを考える。
n=5
k=10^4
p0=runif(k)
bf= (n+1)*p^n
pri=p0/(1-p0)
pos=pri*bf
post_prob=pos/(1+pos)
quantile(post_prob,prob=c(0.025,0.5,0.975))
362:卵の名無しさん
17/12/31 16:32:20.16 dbAzKAtn.net
dataList=list(
a=a,
b=b,
n=n,
z=z,
p=p
)
modelString='
model{
bf = p^z*(1-p)^(n-z)*exp(loggam(a))*exp(loggam(b))/exp(loggam(a+b))*exp(loggam(a+b+n))/exp(loggam(a+z))/exp(loggam(b+n-z))
pri = p0/(1-p0)
pos = pri*bf
post_prob = pos/(1+pos)
p0 ~ dbeta(a,b)
}
'
writeLines(modelString,'TEMPmodel.txt')
jagsModel=jags.model('TEMPmodel.txt', data=dataList)
codaSamples=coda.samples(jagsModel,n.iter=20000,chains=4,var=c('post_prob'))
js=as.matrix(codaSamples)
xlim=c(0,min(1,mean(js[,'post_prob'])*6))
BEST::plotPost(js[,'post_prob'],xlim=xlim,xlab=paste(n,'guess',z,'correct'),cenTend='mean',
showCurve = FALSE,col='gray')
363:卵の名無しさん
17/12/31 17:30:53.02 dbAzKAtn.net
# X[n] ~ N(μ,σ)
# prio.μ ~ N(η,τ)
# post.μ ~ N( [(n/σ2)x_ + (1/τ2)η] /[(n/σ2) + (1/τ2)],
# (σ2/n)*τ2/[(σ2/n) + τ2]
post_norm <- function(eta,tau,n,x_,sigma){
mean=(n/sigma^2*x_ + 1/tau^2*eta)/(n/sigma^2 + 1/tau^2)
sd=(sigma^2/n)*tau^2 / ((sigma^2/n) + tau^2)
return(c(mean,sd))
}
post_norm(eta=180,tau=15,n=5,x_=195,sigma=10)
364:卵の名無しさん
17/12/31 17:31:27.73 dbAzKAtn.net
n=5
X_=195
sigma=10
x=rnorm(n,X_,sigma)
X=(x-mean(x))/sd(x)*sigma + X_
mean(X);sd(X)
eta=180
tau=15
dataList=list(n=n,sigma=sigma,X=X,eta=eta,tau=tau)
modelString='
model{
for(i in 1:n){
X[i] ~ dnorm(mu,1/sigma^2)
}
mu ~ dnorm(eta,1/tau^2)
}
'
writeLines(modelString,'TEMPmodel.txt')
jagsModel=jags.model('TEMPmodel.txt',data=dataList)
codaSamples=coda.samples(jagsModel,n.iter = 100000,var=c('mu'))
js=as.matrix(codaSamples)
hist(js)
mean(js)
var(js)
post_norm(eta=180,tau=15,n=5,x_=195,sigma=10)
365:卵の名無しさん
18/01/01 21:30:04.10 9r6fKQha.net
URLﾘﾝｸ(i.imgur.com)
366:卵の名無しさん
18/01/03 20:47:49.71 CFofwzsi.net
■3囚人問題（英: Thre
367:e Prisoners problem） ある監獄にA、B、Cという3人の囚人がいます 3人のうちランダムに選ばれた1人に恩赦が出ます 誰が恩赦になるかは看守は答えない 囚人Aに看守が「Bは死刑になる」と教えてくれます この時、看守は嘘は言いません 囚人Aに恩赦が与えられる確率は何％でしょうか？
368:卵の名無しさん
18/01/03 20:49:01.81 CFofwzsi.net
死刑囚A,B,CでAが看守に尋ねてBは死刑執行されると告げられたと設定。
恩赦(onsha)を受けるをo,死刑執行されると告(tsuge)げられるをtで表す。
Aが恩赦を受ける確率P(A=o)=1/3
Bが恩赦を受ける確率P(B=o)=1/3
Cが恩赦を受ける確率P(C=o)=1/3
求めたいのは、Bが死刑執行されると告げられた後のAが恩赦を受ける確率P(A=o|B=t)である。
ベイズの公式により
P(A=o|B=t) = P(B=t|A=o)*P(A=o) / P(B=t)
P(B=t) = P(B=t|A=o)*P(A=o) + P(B=t|B=o)*P(B=o) + P(B=t|C=o)*P(C=o)
P(B=t|B=o)=0 Bが恩赦を受けるときBが死刑執行されると告げられる確率=0
P(B=t|C=o)=1 CBが恩赦を受けるときBが死刑執行されると告げられる確率=1
問題は P(B=t|A=o)
恩赦を受けるのがAであるときに看守がCではなくBが死刑執行されると告げる確率は示されていない。
この確率をｐとすると
P(A=o|B=t)は p/(p+1)となる。
もちろんp=0.5であれば、P(A=o|B=t)=1/3と看守に告げられる前と同じである。
ここでｐが一様分布からさまざなβ分布に従うとするとどうなるか、グラフにしてみた。
URLﾘﾝｸ(i.imgur.com)
左の緑が看守がBとCが死刑執行予定であるときにBを選んで答える確率分布。
右の青が看守がBと告げたときのAが恩赦を受ける確率の分布。
369:卵の名無しさん
18/01/03 21:28:14.93 CFofwzsi.net
（タイプミス修正）
死刑囚A,B,CでAが看守に尋ねてBは死刑執行されると告げられたと設定。
恩赦(onsha)を受けるをo,死刑執行されると告(tsuge)げられるをtで表す。
Aが恩赦を受ける確率P(A=o)=1/3
Bが恩赦を受ける確率P(B=o)=1/3
Cが恩赦を受ける確率P(C=o)=1/3
求めたいのは、Bが死刑執行されると告げられた後のAが恩赦を受ける確率P(A=o|B=t)である。
ベイズの公式により
P(A=o|B=t) = P(B=t|A=o)*P(A=o) / P(B=t)
P(B=t) = P(B=t|A=o)*P(A=o) + P(B=t|B=o)*P(B=o) + P(B=t|C=o)*P(C=o)
P(B=t|B=o)=0 Bが恩赦を受けるときBが死刑執行されると告げられる確率=0
P(B=t|C=o)=1 Cが恩赦を受けるときBが死刑執行されると告げられる確率=1
問題は P(B=t|A=o)
恩赦を受けるのがAであるときに看守がCではなくBが死刑執行されると告げる確率は示されていない。
この確率をｐとすると
P(A=o|B=t)は p/(p+1)となる。
もちろんp=0.5であれば、P(A=o|B=t)=1/3と看守に告げられる前と同じである。
ここでｐが一様分布からさまざなβ分布に従うとするとどうなるか、グラフにしてみた。
URLﾘﾝｸ(i.imgur.com)
左の緑が看守がBとCが死刑執行予定であるときにBを選んで答える確率分布。
右の青が看守がBと告げたときのAが恩赦を受ける確率の分布。
370:卵の名無しさん
18/01/03 21:28:54.72 CFofwzsi.net
無情報分布として一様分布を考えると
Aが恩赦を受ける確率の期待値（平均値）は
> 1-log(2)
[1] 0.3068528
となる。
p/(p+1)を [0,1]で定積分すれば求まる。
無情報分布として一様分布を考えると
Aが恩赦を受ける確率の期待値（平均値）は
> 1-log(2)
[1] 0.3068528
となる。
Cが恩赦を受ける確率の期待値（平均値）は
> log(2)
[1] 0.6931472
当然、Bが恩赦を受ける確率は０ 👀
Rock54: Caution(BBR-MD5:0be15ced7fbdb9fdb4d0ce1929c1b82f)
371:卵の名無しさん
18/01/04 16:17:58.24 nNHCcGvL.net
P(A=o|B=t) = P(B=t|A=o)*P(A=o) / P(B=t)
P(B=t) = P(B=t|A=o)*P(A=o) + P(B=t|B=o)*P(B=o) + P(B=t|C=o)*P(C=o) = (p + q+ 1-q)*1/3 = (p + 1)*1/3
P(B=t|A=o)*P(A=o) = p*1/3
P(B=t|B=o)*P(B=o) = q*1/3
P(B=t|C=o)*P(C=o) = (1-q)*1/3
P(A=o|B=t) = p / (p +1)
372:卵の名無しさん
18/01/05 09:01:44.41 HIwtxZA8.net
> vonNeumann
function(PDF,xmin=0,xmax=1,N=10000,Print=TRUE,...){
xx=seq(xmin,xmax,length=N+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=4)
invisible(Rand)
}
373:卵の名無しさん
18/01/05 13:47:38.14 2i23lFVy.net
####
stanString2='
functions{
real jisaku_lpdf(real y, real a, real b){
real temp;
temp = a*y + b;
return log(temp);
}
}
data{
real a;
real b;
}
parameters{
real<lower=0,upper=1> p;
}
transformed parameters{
real q;
q = p/(p+1);
}
model{
target += jisaku_lpdf( p | a, b);
}
'
374:卵の名無しさん
18/01/06 05:46:40.68 Qus13FVj.net
amnesty <- function(p) p/(p+1)
f=function(ab,k=10^4,ma=1/3){
a=ab[1]
b=ab[2]
rn=rbeta(k,a,b)
(mean(amnesty(rn)) - ma)^2
}
optim (c(1,1),f, method='Nelder')
optim (c(2,2
375:),f, method='Nelder')
376:卵の名無しさん
18/01/07 20:56:03.58 63fon0Wx.net
■3囚人問題（英: Three Prisoners problem）
ある監獄にA、B、Cという3人の囚人がいます
3人のうちランダムに選ばれた1人に恩赦が出ます
誰が死刑になるかは看守は本人には答えない
囚人Aに看守が「Bは死刑になる」と教えた
この時、看守は一定の確率で嘘をつく。
BとCが死刑になるときは50%の確率でBと答える。
囚人Aに恩赦が与えられる確率は何％でしょうか？
377:卵の名無しさん
18/01/09 19:32:17.33 WXFR1iXr.net
■3囚人問題（英: Three Prisoners problem）
ある監獄にA、B、Cという3人の囚人がいます。
3人のうちランダムに選ばれた1人に恩赦が出ます。
誰が死刑になるかは看守は決して本人には教えない。
囚人AがB、Cの少なくともどちらかは死刑になるのだから教えてくれと看守に尋ねた。看守が本人に教えるのではないので「Bは死刑になる」とAに教えた。
この時、看守は1/3の確率で嘘をつく。
BとCが死刑になるときは50%の確率でBと答える。
囚人A、B、Cに恩赦が与えられる確率はそれぞれ何％でしょう？
378:卵の名無しさん
18/01/09 19:41:09.46 WXFR1iXr.net
>>447
p : P(B=t|A=o)Aが恩赦（BとCが死刑執行される）とき看守がBと答える確率
q : 看守が嘘をつく確率
P(B=t|B=o) Bが恩赦を受けるときBが死刑執行されると告げられる確率 = q
P(B=t|C=o) Cが恩赦を受けるときBが死刑執行されると告げられる確率 = 1-q
P(A=o|B=t) = P(B=t|A=o)*P(A=o) / P(B=t)
P(B=t) = P(B=t|A=o)*P(A=o) + P(B=t|B=o)*P(B=o) + P(B=t|C=o)*P(C=o)
= p * P(A=o) + q * P(B=o) + (1-q) * P(C=o)
P(A=o|B=t) = p*P(A=o) / ( p*P(A=o) + q * P(B=o) + (1-q) * P(C=o) )
P(A=o)= P(B=o)= P(C=o) = 1/3ならば
P(A=o|B=t) = p /(p+1)
P(B=o|B=t) = q/(p+1)
P(C=o|B=t) = (1-q)/(p+1)
379:卵の名無しさん
18/01/18 15:55:28.61 mFfO4JsF.net
プレーヤーが選んだ箱をA、司会者モンティーホールが開けたハズレの箱をB、残った箱をCとする。
Aがアタリ(atari)の確率をP(A=a)
司会者がBを開ける(open)確率をP(B=o)と表すことにする。
残った箱Cがアタリである確率P(C=a|B=o)は
ベイズの公式から
P(C=a|B=o) = (P(B=o|C=a)P(C=a)) / P(B=o)
P(B=o) = P(B=o|A=a)P(A=a) + P(B=o|B=a)P(B=a) + P(B=o|C=a)P(C=a)
= P(B=o|A=a)*1/3 + 0*1/3 + 1*1/3 ここで P(B=o|A=a)=pとおくと
= p*1/3 + 1/3
ゆえに
P(C=a|B=o) = (P(B=o|C=a)P(c=a)) / P(B=o) = (1*1/3) / (p*1/3 + 1/3) = 1/(p+1)
380:卵の名無しさん
18/01/18 15:58:49.78 mFfO4JsF.net
BがハズレとわかったあとでAがアタリである確率
P(A=a|B=o) = P(B=o|A=a)P(A=a)/P(B=o)
P(B=o) = P(B=o|A=a)P(A=a) + P(B=o|B=a)P(B=a) + P(B=o|C=a)P(C=a)
P(B=o|A=a)はAがアタリであるときにBがハズレとして開けられる確率pは問題で示されていない。
不十分理由の原則に準じてpを０．５とするか一様分布に従うとするのが一般的だと思う。
P(B=o|A=a)=pとおくと
P(B=o) = P(B=o|A=a)P(A=a) + P(B=o|B=a)P(B=a) + P(B=o|C=a)P(C=a)
= p*1/3 + 0*1/3 + 1*1/3
= p*1/3 +1/3
ゆえに
P(A=a|B=o) = (p*1/3) / ( p*1/3 + 1/3 ) = p/(p+1)となる。
p=0.5ならBがハズレというデータはAがあたりの確率に影響を与えず1/3である。
381:卵の名無しさん
18/01/18 15:59:07.36 mFfO4JsF.net
p: Aがアタリの時に司会者がBを開ける確率
P(A=a|B=o) = p/(p+1) Bが開けられた後、Aがアタリの確率 (1)
P(C=a|B=o) = 1/(p+1) Bが開けられた後、Cがアタリの確率 (2)
(1)/(2) = p なので (2)は(1)以上である。(∵0<= p <=1)
ゆえに
残った箱Cの方がアタリの確率は高い。
382:卵の名無しさん
18/01/21 07:14:29.94 VjwevIsc.net
ゴルゴ15は1発1中
とする。
各々10000発撃ったときゴルゴ15の命中数の期待値はいくらか？
確率密度とかベータ分布とかを使わずに説明するなら、重み付き平均という考え方で説明するしかないかな？
命中率が0.5なら2回に1回は1発1中（確率0.5)
命中率が0.8なら10回に8回は1発1中（確率0.8)
となる。
体重100ｋｇの牛が100頭
体重99kgの牛が99頭
体重98kgの牛が98頭
・・・
体重2kgの牛が2頭
体重1kgの牛が1頭
牛の平均体重の計算と同じ
n=100
x=seq(0,1,length=n+1)
sum(x*x/sum(x))
sum(x^2)/sum(x)
2/3
(sum_x=n*(n+1)/2/n) # (n+1)/2
(sum_x2=n*(n+1)*(2*n+1)/6/(n^2)) # (n+1)*(2*n+1)/n/6
sum_x2/sum_x # (2*n+1)/n/3 = 2/3+1/3/n
n→∞
で2/3に集束する。 命中数の期待値は10000*2/3=6667
383:卵の名無しさん
18/03/08 13:29:28.47 1C4m5UpT.net
URLﾘﾝｸ(i.imgur.com)
リリースされているデータだと効き目はこんなもん
384:卵の名無しさん
18/03/08 17:51:23.53 jJMnC74H.net
R でログランク検定を実行するには，観察時間を示す変数を time，打ち�
385:ﾘりフラグを event，グループを group とすれば，survdiff(Surv(time,event)~group) とすればよい。この例の場合なら，下枠内の通り。 なお，一般化ウィルコクソン検定をするには survdiff(Surv(time,event)~group,rho=1) とすればよい。 it13-4-2006.R require(survival) time <- c(4,6,8,9,5,7,12,14) event <- c(1,1,1,1,1,1,1,1) group <- c(1,1,1,1,2,2,2,2) survdiff(Surv(time,event)~group)
386:卵の名無しさん
18/03/08 21:28:51.57 jJMnC74H.net
サイコロを３０００回振って１の目が４９０回でたサイコロはイカサマサイコロか？
１０％までの歪みは許容する。
387:卵の名無しさん
18/03/08 21:29:24.44 jJMnC74H.net
crooked <- function(n,r,H0=0.5,d=0.1,credMass=0.95,Print=TRUE,...){
hdi=HDInterval::hdi(qbeta,credMass,shape1=1+r,shape2=1+n-r)
return(hdi)
}
crooked(3000,490,H0=1/6,xlim=c(0.12,0.20))
binom::binom.confint(490,3000)
1/6*c(0.9,1.1)
グラフ化するとHDI ⊂ ROPEになっている。
URLﾘﾝｸ(i.imgur.com)
388:卵の名無しさん
18/03/12 07:04:43.41 RtLIiZVI.net
# URLﾘﾝｸ(www.ncbi.nlm.nih.gov)
test\disease present absent
pos TP(a) FN(b)
neg FP(c) TN(d)
a+c=33
b+d=15
a+b=33
a/(a+c)=0.697
d/(b+d)=0.333
a/33=0.697
a=23
d/15=0.333
d=5
c=10
b=10
test\disease present absent
pos 23 TP(a) 10 FN(b)
neg 10 FP(c) 5 TN(d)
# URLﾘﾝｸ(statpages.info)
389:卵の名無しさん
18/03/12 13:24:27.62 UNJR7sdw.net
頭がある という所見は髄膜炎の診断に感度１００％である。
角がある という所見は髄膜炎の診断に特異度１００％である。
俺の経験上、こういう話をシリツ医大卒に振っても興味を示す奴はいないね。
390:卵の名無しさん
18/03/20 23:04:27.14 ql3rAyFF.net
pLR=.99/.01
nLR=.01/.99
preo=0.001/(1-0.001)
preo=0.1/(1-0.1)
poso=preo*pLR
poso/(1+poso)
391:卵の名無しさん
18/03/21 20:43:38.90 5ZNg8kVw.net
# modified Wald by Agresti and Coull
mWald <- function(S,n,cl=0.95){
z=qnorm(1-(1-cl)/2)
n_tilde=n+z^2
p_tilde=(S+z^2/2)/n_tilde
W=z*sqrt(p_tilde*(1-p_tilde)/n_tilde)
CI=c(mean=S/n,lower=p_tilde-W,upper=p_tilde+W)
return(CI)
}
binom::binom.agresti.coull(S,n)
392:卵の名無しさん
18/03/22 08:58:41.40 IaVjTmR+.net
3/n ≒ 1 - n√0.05
393:卵の名無しさん
18/03/23 10:37:59.90 OxTHnFDc.net
#標準偏差の信頼区間
sdCI=function(x,conf.level=0.95){
n=length(x)
df=n-1
sd=sd(x)
lwr=qchisq((1-conf.level)/2,df)
upr=qchisq((1-conf.level)/2,df,lower.tail=FALSE)
CI=data.frame(lower=sqrt((n-1)/upr)*sd,sd=sd,upper=sqrt((n-1)/lwr)*sd )
return(CI)
}
##差の信頼区間(生データなし)
DifCI=function(n1,n2,m1,m2,sd1,sd2){
pooledV=((n1-1)*sd1^2+(n2-1)*sd2^2)/(n1-1+n2-1)
SE12=sqrt((1/n1+1/n2)*pooledV)
w=qt(.975,n1-1+n2-1)*SE12
ci=c(m1-m2-w,m1-m2+w)
names(ci)=c("lower","upper")
return(ci)
}
# 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
pt(abs(T),n1-1+n2-1,lower.tail = FALSE)
}
394:卵の名無しさん
18/03/23 10:38:50.82 OxTHnFDc.net
## m sd SEM n
non=c(0.52,0.25,0.027,88)
rec=c(0.38,0.32,0.034,89)
eli=c(0.40,0.26,0.048,28)
# 生データなしで分散分析
lh=rbind(non,rec,eli)
colnames(lh)=c("m","sd","SEM","n") ; lh
mean.G=sum(lh[,"m"]*lh[,"n"])/sum(lh[,"n"])
SS.bit=sum((lh[,"m"]-mean.G)^2*lh[,"n"])
SS.wit=sum(lh[,"sd"]^2*(lh[,"n"]-1))
df.bit=nrow(lh)-1
df.wit=sum(lh[,"n"]-1)
MS.bit=SS.bit/df.bit
MS.wit=SS.wit/df.wit
F.ratio=MS.bit/MS.wit
pf(F.ratio,df.bit,df.wit,lower.tail=FALSE) # 0.003720507
(η2=(SS.bit)/(SS.bit+SS.wit)) # 0.05387927
395:卵の名無しさん
18/03/23 10:40:13.81 OxTHnFDc.net
#prop.test(S,n)$conf.int # S:success, n:number of trial
prop=function(S,n)prop.test(S,n)$conf.int[1:2]
#binom.test(S,n)$conf.int # S:success, n:number of trial
binom=function(S,n) binom.test(S,n)$conf.int[1:2]
## S:success, n:number of trial
binomCI=function(S,n,conf.level=0.95){
upper=uniroot(f=function(x)pbinom(S,n,x)-(1-conf.level)/2,c(0,1))$root
lower=uniroot(f=function(x)pbinom(S,n,x,lower.tail=FALSE)-(1-conf.level)/2,c(0,1))$root
CI=data.frame(lower=lower, upper=upper)
return(round(CI,3))}
##
wCI=function(S,n,cl=0.95){ #modified Wald CI for Proportion
z=qnorm(1-(1-cl)/2)
p.=(S+0.5*z^2)/(n+z^2)
W=z*sqrt((p.*(1-p.))/(n+z^2))
C.I.=data.frame(lower=p.-W,upper=p.+W)
return(C.I.)
}
##
bCI<-function(r,n,cl=0.95){
p<-r/n
q<-1-p
sd<-sqrt (n*p*q)
SE<-sd/sqrt(n)
Z<-qnorm(1-(1-cl)/2)
ci<-data.frame(lower=p-Z*SE/sqrt(n),upper=p+Z*SE/sqrt(n))
return(ci)
}
396:卵の名無しさん
18/03/25 14:44:59.55 Z8Ybstvj.net
Signif?cosis. Manifested by a failure to discern between
biological and statistical signif?cance (6). Individuals with
signif?cosis fail to realize that just because something is
unlikely to have occurred by chance doesn’t mean it’s
important (17). See also Borderline Probability Disorder.
Borderline Probability Disorder. Afflicted individuals
may dismiss the potential importance of results with P=0.06
while unquestioningly accepting the importance of results with P=0.05
URLﾘﾝｸ(journals.asm.org)
397:卵の名無しさん
18/03/25 16:56:05.47 Z8Ybstvj.net
r1=10 ; r2=16
n1=5936 ; n2=3403
a=r1
b=n1-r1
c=r2
d=n2-r2
Fragile <- function(a,b,c,d){
mat=matrix(c(a,b,c,d),2,2,byrow=TRUE,dimnames=list(c("Treated","Control"),c("Event","No Event")))
print(mat)
p=fisher.test(mat)$p.value
cat(paste('p-value =',round(p,5)))
i=0
while(p < 0.05 ){
p=fisher.test(matrix(c(a+i,b-i,c,d),2,2,byrow=TRUE))$p.value
i=i+1
}
cat('\n','\n')
cat(paste('Fragile Index = ',i,'\n','\n'))
print(matrix(c(a+i,b-i,c,d),2,2,byrow=TRUE,dimnames=list(c("Treated","Control"),c("Event","No Event"))))
cat(paste('p-value =',round(p,5)))
invisible(i)
}
398:卵の名無しさん
18/03/26 07:18:15.20 HyN4R5wr.net
# High Normal Total
# Event a b x
# No Event c d y
sqrt(1.37*4.09)
x+y=43407
x=1280
y=43407-1280
HR=x/(x+y) ; HR
HR1= 2.37*HR ; HR1
a/(a+c)=HR1
c=(1-HR1)*a/HR1
V(logRR)=1/a-1/x+1/c-1/y
SE(logRR)=sqrt(1/a-1/x+1/c-1/y)
HR1*exp(1.96*SE(logRR))=4.09
399:卵の名無しさん
18/03/26 07:19:05.61 HyN4R5wr.net
library(fmsb)
hazardratio <-function (a, b, PT1, PT0, conf.level = 0.95) {
.M <- a + b
.T <- PT1 + PT0
.MAT <- matrix(c(a, b, .M, PT1, PT0, .T), 3, 2)
colnames(.MAT) <- c("Cases", "Person-time")
rownames(.MAT) <- c("Exposed", "Unexposed", "Total")
class(.MAT) <- "table"
# print(.MAT)
ESTIMATE <- (a/PT1)/(b/PT0)
norm.pp <- qnorm(1 - (1 - conf.level)/2)
.CHI <- (a - (PT1/.T) * .M)/sqrt(.M * (PT1/.T) * (PT0/.T))
p.v <- 2 * (1 - pnorm(abs(.CHI)))
RRL <- ESTIMATE * exp(-norm.pp * sqrt(1/a + 1/b))
RRU <- ESTIMATE * exp(norm.pp * sqrt(1/a + 1/b))
CINT <- c(RRL, RRU)
attr(CINT, "conf.level") <- conf.level
RVAL <- list(p.value = p.v, conf.int = CINT, estimate = ESTIMATE,
method = "Incidence rate ratio estimate and its significance probability",
data.name = paste(deparse(substitute(a)), deparse(substitute(b)),
deparse(substitute(PT1)), deparse(substitute(PT0))))
class(RVAL) <- "htest"
return(RVAL)
}
f <- function(a,u0=1.37) hazardratio(a, x, a/HR1, x+y)$conf[1] - u0
uniroot(f,c(1,1000))$root
a=13
(1-HR1)*a/HR1
c=173
400:卵の名無しさん
18/03/26 16:51:22.07 xQgYtict.net
Run any R code you like. There are over three thousand R packages preloaded.
URLﾘﾝｸ(rdrr.io)
401:卵の名無しさん
18/03/27 06:27:14.34 7KQq5UQg.net
N=1000
X=rnorm(N,0,1)
Y=rnorm(N,1,1)
n=17
f <- function(){
x=sample(X,n)
y=sample(Y,n)
x-y
}
d=replicate (10^5, f())
hist(d)
g <- function(){
x=sample(X,n)
y=sample(Y,n)
t.test(x,y)$p.value
}
p=replicate (10^4, g())
hist(p)
mean(p<0.05)
402:卵の名無しさん
18/03/27 06:36:24.73 7KQq5UQg.net
N=1000
X=rnorm(N,0,1)
Y=rnorm(N,1,1)
n=16
f <- function(){
x=sample(X,n)
y=sample(Y,n)
x-y
}
d=replicate (10^5, f())
hist(d)
g <- function(){
x=sample(X,n)
y=sample(Y,n)
t.test(x,y)$p.value
}
p=replicate (10^4, g())
hist(p)
hist(-log10(p))
mean(log10(p)<log10(0.05))
403:卵の名無しさん
18/03/27 06:42:13.61 7KQq5UQg.net
N=1000
X=rnorm(N,0,1)
Y=rnorm(N,1,1)
n=16
f <- function(){
x=sample(X,n)
y=sample(Y,n)
mean(x)-mean(y)
}
d=replicate (10^5, f())
hist(d)
sd(d)
g <- function(){
x=sample(X,n)
y=sample(Y,n)
t.test(x,y)$p.value
}
p=replicate (10^5, g())
hist(p)
hist(-log10(p))
mean(log10(p)<log10(0.05))
404:卵の名無しさん
18/03/27 06:42:36.33 7KQq5UQg.net
URLﾘﾝｸ(rsos.royalsocietypublishing.org)
405:卵の名無しさん
18/03/27 06:50:47.55 7KQq5UQg.net
> N=1000
> X=rnorm(N,0,1)
> Y=rnorm(N,1,1)
> n=16
> f <- function(){
+ x=sample(X,n)
+ y=sample(Y,n)
+ mean(x)-mean(y)
+ }
> d=replicate (10^5, f())
> hist(d)
> sd(d)
[1] 0.3574393
>
> g <- function(){
+ x=sample(X,n)
+ y=sample(Y,n)
+ t.test(x,y)$p.value
+ }
> p=replicate (10^5, g())
> hist(p)
> hist(-log10(p))
> mean(log10(p)<log10(0.05))
[1] 0.79713
>
>
406:卵の名無しさん
18/03/27 06:52:28.91 7KQq5UQg.net
> power.t.test(n=16,delta=1)
Two-sample t test power calculation
n = 16
delta = 1
sd = 1
sig.level = 0.05
power = 0.7813965
alternative = two.sided
NOTE: n is number in *each* group
407:卵の名無しさん
18/03/27 08:07:28.48 7KQq5UQg.net
N=1000
X=rnorm(N,0,1)
n=16
f <- function(){
x=sample(X,n)
y=sample(X,n)
mean(x)-mean(y)
}
d=replicate (10^5, f())
hist(d)
sd(d)
g <- function(){
x=sample(X,n)
y=sample(X,n)
t.test(x,y)$p.value
}
p=replicate (10^4, g())
hist(p)
hist(-log10(p))
mean(log10(p)<log10(0.05))
408:卵の名無しさん
18/03/27 10:32:56.09 7KQq5UQg.net
>>385
理論上、
mean(p<0.05)は
power.t.test(n=16,delta=1)$powerに一致。
409:卵の名無しさん
18/03/28 09:07:02.26 lnrrtf8h.net
calc.FPR = function(nsamp,pval,sigma,prior,delta1){
# sdiff=sqrt(sigma^2/nsamp + sigma^2/nsamp)
ns1=nsamp
ns2=nsamp
# CALC sdiff = sd of difference between means
sdiff=sqrt(sigma^2/ns1 + sigma^2/ns2)
df=ns1+ns2-2
# Note FPR doesn't need calculation of power for p-equals case
#
#under H0, use central t distribution
tcrit=qt((1-pval/2),df,ncp=0)
x0=tcrit
y0=dt(x0,df,0)
#
# under H1 use non-central t distribution
ncp1=delta1/sdiff #non-centrality paramater
x1=x0 #tcrit
y1=dt(x1,df,ncp=ncp1)
# check solution
# pcheck=pt(y1,df,ncp=ncp1)
# pcheck
# Calc false positive risk
p0=2*y0
p1=y1
FPR=((1-prior)*p0)/(((1-prior)*p0) + prior*p1)
FPR
output=c(FPR,x0,y0,x1,y1)
return(output)
}
# end of function calc.FPR
#
410:卵の名無しさん
18/03/28 09:07:17.56 lnrrtf8h.net
# calc.FPR0 gives FPR for given nsamp, for p-less-than case
calc.FPR0 = function(nsamp,pval,sigma,prior,delta1){
myp=power.t.test(n=nsamp,sd=sigma,delta=delta1,sig.level=pval,type="two.sample",alternative="two.sided",power=NULL)
power = myp$power
PH1=prior
PH0=1-PH1
FPR0=(
411:pval*PH0/(pval*PH0 + PH1*power)) output=c(FPR0,power) return(output) } #
412:卵の名無しさん
18/03/28 09:27:44.91 lnrrtf8h.net
# False Positive Report Probability
FPRP <- function(n, p.value, prior, effect.size=1){
power=power.t.test(n=n,delta=effect.size,sig.level=p.value)$power
FP=(1-prior)*p.value
TP=prior*power
return(FP/(TP+FP))
}
FPRP(16,p.value=0.045,prior=0.10)
# サンプルサイズが違うとき
FPRP2 <- function(n1,n2, p.value, prior, effect.size=1,...){
power=pwr::pwr.t2n.test(n1, n2, d=effect.size, sig.leve=p.value, ...)$power
FP=(1-prior)*p.value # p-less-than case
TP=prior*power
return(FP/(TP+FP))
}
FPRP2(10,6,p.value=0.05,prior=0.05)
413:卵の名無しさん
18/03/28 10:44:32.63 lnrrtf8h.net
# False Positive Report Probability
FPRP <- function(n, p.value, prior, effect.size=1){
power=power.t.test(n=n,delta=effect.size,sig.level=p.value)$power
FP=(1-prior)*p.value
TP=prior*power
return(FP/(TP+FP))
}
FPRP(16,p.value=0.045,prior=0.10)
# サンプルサイズが違うとき
FPRP2 <- function(n1,n2, p.value, prior, effect.size=1,...){
power=pwr::pwr.t2n.test(n1, n2, d=effect.size, sig.leve=p.value, ...)$power
FP=(1-prior)*p.value # p-less-than case
TP=prior*power
return(FP/(TP+FP))
}
FPRP2(10,6,p.value=0.05,prior=0.05)
414:卵の名無しさん
18/03/28 18:20:29.83 izUkh1ln.net
更新版の
URLﾘﾝｸ(rsos.royalsocietypublishing.org)
を読んでたらRのスクリプトが公開されていた。
URLﾘﾝｸ(figshare.com)
415:卵の名無しさん
18/03/28 20:36:10.95 izUkh1ln.net
FPRP = P(H0|Rejected) = P(Rejected|H0)P(H0) / [ P(Rejected|H0)P(H0) + P(Rejected|H1)P(H1)]
= alpha*(1-prior) / [ alpha*(1-prior) + power*prior ] (1)
= p.value*(1-prior) / [ p.value*(1-prior) + power*prior ] (2)
416:卵の名無しさん
18/03/29 14:32:00.14 HGAAnleD.net
シミュレーションデータ
> length(CONTROLS) ; length(HYPERTENSION)
[1] 17
[1] 18
> mean(CONTROLS) ; mean(HYPERTENSION)
[1] 263
[1] 257
> sd(CONTROLS) ; sd(HYPERTENSION)
[1] 87
[1] 59
で条件に合致。
> t.test(CONTROLS,HYPERTENSION,var=TRUE)
Two Sample t-test
data: CONTROLS and HYPERTENSION
t = 0.24003, df = 33, p-value = 0.8118
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-44.85718 56.85718
JAGSを使ってMCMCしてみた(パッケージBESTを改造して使用）。
URLﾘﾝｸ(i.imgur.com)
理論値より信頼区間幅が広くなるが、まあ、シミュレーションなのでこんなもんだろう。
> HDIofMCMC(muDiff)
[1] -48.86974 61.99241
自分で計算してグラフ化しながら読み進むのは楽しい。
417:卵の名無しさん
18/03/29 22:09:13.07 HGAAnleD.net
library(pwr)
args(pwr.t2n.test)
n1=17;n2=18;m1=263;m2=257;sd1=87;sd2=59
sp_sq = ((n1-1)*sd1^2 + (n2-1)*sd2^2)/(n1-1+n2-1) # pooled variance
sp=sqrt(sp_sq) ; sp # pooled sd
sep=sqrt(1/n1+1/n2)*sp # se for pooled sd
f20.1 <- function(x) pwr.t2n.test(n1=17,n2=18, d=x, sig.level = 0.05)$power
xx=0:150
plot(xx,sapply(xx/sp,f20.1),type='l',ylab='power',xlab='mean difference')
f20.1(0)
f20.1(1)
#
z1=383;z2=373;n1=7685;n2=7596
p1=z1/n1
p2=z2/n2
p1/p2
ES.h(p1,p2)
g20.1 <- function(x){
pwr.2p2n.test(h=ES.h(p1,x*p1),n1=n1,n2=n2,sig.level=0.05)$power
}
g20.1(1)
g20.1(0.7)
rr=seq(0,1,by=0.01)
plot(rr,sapply(rr,g20.1),type='l',xlab='relative risk ratio',ylab='power')
418:卵の名無しさん
18/03/30 06:00:08.04 oNlAm/3j.net
# False Positive Report Probability for Proportion
FPRPP <- function(r1,r2,n1,n2, prior){
p1=r1/n1
p2=r2/n2
p.value=prop.test(c(r1,r2),c(n-1,n2))$p.value
pwr::pwr.2p2n.test(h=ES.h(p1,x*p1),n1=n1,n2=n2,sig.level=p.value)$power
FP=(1-prior)*p.value # p-less-than case
TP=prior*power
return(FP/(TP+FP))
}
r1=10 ; r2=16
n1=5936 ; n2=3403
FPRPP(10,16,5936,3403,0.5)
419:卵の名無しさん
18/03/30 06:14:57.20 oNlAm/3j.net
# False Positive Report Probability for Proportion
FPRPP <- function(r1,r2,n1,n2, prior){
p.value=prop.test(c(r1,r2),c(n1,n2))$p.value
p1=r1/n1
p2=r2/n2
h=pwr::ES.h(p1,p2)
power=pwr::pwr.2p2n.test(h=h,n1=n1,n2=n2,sig.level=p.value)$power
FP=(1-prior)*p.value # p-less-than case
TP=prior*power
return(FP/(TP+FP))
}
r1=10 ; r2=16
n1=5936 ; n2=3403
FPRPP(r1,r2,n1,n2, prior=0.5)
420:卵の名無しさん
18/03/30 07:49:47.13 oNlAm/3j.net
URLﾘﾝｸ(clincalc.com)
421:卵の名無しさん
18/03/30 09:44:46.63 oNlAm/3j.net
# Fragile Index
FragileIndex<- function(r1,r2,n1,n2){
if(
422:r1/n1 < r2/n2){ a=r1 b=n1-r1 c=r2 d=n2-r2 }else{ a=r2 b=n2-r2 c=r1 d=n1-r1 } mat=matrix(c(a,b,c,d),2,2,byrow=TRUE,dimnames=list(c("low","high"),c("Event","No Event"))) print(mat) p=fisher.test(mat)$p.value cat(paste('p-value =',round(p,5))) i=0 while(p < 0.05 ){ p=fisher.test(matrix(c(a+i,b-i,c,d),2,2,byrow=TRUE))$p.value i=i+1 } cat('\n','\n') cat(paste('Fragile Index = ',i,'\n','\n')) print(matrix(c(a+i,b-i,c,d),2,2,byrow=TRUE,dimnames=list(c("low+FI","high-FI"),c("Event","No Event")))) cat(paste('p-value =',round(p,5))) invisible(i) }
423:卵の名無しさん
18/03/30 10:30:37.99 oNlAm/3j.net
# URLﾘﾝｸ(www.jclinepi.com)(13)00466-6/pdf
# Fragile Index
FragileIndex<- function(r1,r2,n1,n2,Fisher=TRUE){
if(r1/n1 < r2/n2){
a=r1
b=n1-r1
c=r2
d=n2-r2
}else{
a=r2
b=n2-r2
c=r1
d=n1-r1
}
mat=matrix(c(a,b,c,d),2,2,byrow=TRUE,dimnames=list(c("low","high"),c("Event","No Event")))
print(mat)
FUN=ifelse(Fisher,fisher.test,chisq.test)
p=FUN(mat)$p.value
cat(paste('p-value =',round(p,5)))
i=0
while(p < 0.05 ){
p=fisher.test(matrix(c(a+i,b-i,c,d),2,2,byrow=TRUE))$p.value
i=i+1
}
cat('\n','\n')
cat(paste('Fragile Index = ',i,'\n','\n'))
print(matrix(c(a+i,b-i,c,d),2,2,byrow=TRUE,dimnames=list(c("low+FI","high-FI"),c("Event","No Event"))))
cat(paste('p-value =',round(p,5)))
invisible(i)
}
424:卵の名無しさん
18/03/30 14:34:22.51 oNlAm/3j.net
rm(list=ls(all=TRUE))
graphics.off()
source('URLﾘﾝｸ(raw.githubusercontent.com)')
BMV_vs_ETI <- function(zi=44 , Ni=1018 , zj=43 , Nj=1022, ROPE=c(-0.005,0.005),COLOR='skyblue'){
y=c(rep(1,zi),rep(0,Ni-zi),rep(1,zj),rep(0,Nj-zj))
s=rep(1:2,c(Ni,Nj))
shape1=1 ; shape2=1 # JAGS prior : beta(shape1,shape2)
myData=data.frame(y=y,s=s)
Ntotal = length(y)
Nsubj = length(unique(s))
dataList = list(
y = y ,
s = s ,
Ntotal = Ntotal ,
Nsubj = Nsubj
)
# JAGS model
modelString = paste0("
model {
for ( i in 1:Ntotal ) {
y[i] ~ dbern( theta[s[i]] )
}
for ( sIdx in 1:Nsubj ) {
theta[sIdx] ~ dbeta(", shape1,',' , shape2," )
}
}
")
# end of modelString
425:卵の名無しさん
18/03/30 14:34:54.83 oNlAm/3j.net
writeLines( modelString , con="TEMPmodel.txt" )
jagsModel = jags.model( file="TEMPmodel.txt" , data=dataList)
update( jagsModel)
codaSamples = coda.samples( jagsModel , variable="theta", n.iter=10000 )
mcmcMat=as.matrix(codaSamples)
dif=mcmcMat[,1]-mcmcMat[,2]
plotPost(dif,ROPE=ROPE,compVal=0,xlab=quote(Delta),breaks=30,cenT='mean',col=COLOR)
invisible(dif)
}
# favorable functional survival at day 28
BMV_vs_ETI(zi=44 , Ni=1018 , zj=43, Nj=1022, ROPE=c(-0.01,0))
# Survival to hospital admission 294/1018 vs 333/1022
BMV_vs_ETI(294,1018,333,1022,ROPE=c(-0.01,0.01))
# global survival at day 28 55/1018 vs 54/1022
BMV_vs_ETI(55,1018,54,1022,ROPE=c(-0.01,0.01))
426:卵の名無しさん
18/03/30 22:20:19.00 /bcqdZEY.net
Generic=rlnorm(10,4.5,0.2)
Brand=rlnorm(20,4.3,0.1)
while(!Equiv(Generic,Brand)){
Generic=rlnorm(10,4.5,0.2)
Brand=rlnorm(20,4.3,0.1)
}
Generic
Brand
t.test(Generic,Brand)$p.value
tt=t.test(Generic,conf.level = 0.90) ; tt
lo=tt$conf[1]
hi=tt$conf[2]
m=mean(Brand)
0.80*m < lo
hi < 1.25*m
427:卵の名無しさん
18/03/31 07:28:51.22 7NaQ5iyd.net
# How big the sample should be to show non-inferiority
r1=44; n1=1018;r2=43;n2=1022
BMV_vs_ETI(r1,n1,r2,n2)
risk.difference(r1,r2,n1,n2)
f <- function(x,margin=0.01){
rd=fmsb::riskdifference(x*r1,x*r2,x*n1,x*n2)
lo=rd$conf[1]
hi=rd$conf[2]
- margin - lo
}
m=uniroot(f,c(0.001,5))$root ; m
BMV_vs_ETI(round(m*r1),round(m*n1),round(m*r2),round(m*n2),COLOR = 'pink')
428:卵の名無しさん
18/04/01 20:44:41.26 yXplTzVs.net
Freedman50 <- function(){
X=matrix(rnorm(50*100),ncol=100) # X[50 variables,100 samples]
Y=rnorm(100)
ans=lm(Y~ X[1,]+X[2,]+X[3,]+X[4,]+X[5,]+X[6,]+X[7,]+X[8,]+X[9,]+X[10,]+X[11,]+X[12,]+X[13,]+X[14,]+X[15,]+X[16,]+X[17,]+X[18,]+X[19,]+X[20,]+
X[21,]+X[22,]+X[23,]+X[24,]+X[25,]+X[26,]+X[27,]+X[28,]+X[29,]+X[30,]+X[31,]+X[32,]+X[33,]+X[34,]+X[35,]+X[36,]+X[37,]+X[38,]+X[39,]+X[40,]+
X[41,]+X[42,]+X[43,]+X[44,]+X[45,]+X[46,]+X[47,]+X[48,]+X[49,]+X[50,])
anova.ans=anova(ans)
p.value=anova.ans$`Pr(>F)`[1] # p.value
429: of regression smry=summary(ans) r.squared=smry$r.squared pv.coef=smry$coef[-1,4] pv.coef025=sum(pv.coef<0.25) # p.value of coefficient pv.coef005=sum(pv.coef<0.05) data.frame(p.value, r.squared,pv.coef025,pv.coef005) } Freedman50() res50=replicate(10^3,Freedman50()$p.value) BEST::plotPost(res50,compVal = 0.05)
430:卵の名無しさん
18/04/02 09:48:50.18 qPQfofGo.net
f<- function(d,beta=.80,alpha=.05) 2*(abs((qnorm(1-alpha/2))+abs(qnorm(1-beta)))/d)^2
f(5/10)
power.t.test(delta=0.5,power=0.80)$n
dd=seq(0.01,0.99,by=0.01)
g<- function(x) f(x,.50)/f(x,.80)
sapply(dd,g)
(2*(abs(1.96+abs(qnorm(1-0.80))))^2)/(2*(abs(1.96+abs(qnorm(1-0.50))))^2)
(1.96+0.84)^2/(1.96+0)^2
431:卵の名無しさん
18/04/02 10:00:29.14 qPQfofGo.net
f<- function(d,beta=.80,alpha=.05) 2*(abs((qnorm(1-alpha/2))+abs(qnorm(1-beta)))/d)^2
f(5/10)
power.t.test(delta=0.5,power=0.80)$n
dd=seq(0.01,0.99,by=0.01)
g<- function(x) f(x,.50)/f(x,.80)
sapply(dd,g)
(2*(abs(1.96+abs(qnorm(1-0.80))))^2)/(2*(abs(1.96+abs(qnorm(1-0.50))))^2)
(1.96+0.84)^2/(1.96+0)^2
dd=seq(0.01,0.99,by=0.01)
h<-function (x){
power.t.test(d=x, power=0.80)$n/power.t.test(d=x, power=0.5)$n
}
res=sapply (dd,h)
plot(dd,res)
summary (res)
432:卵の名無しさん
18/04/03 00:15:01.95 wYZjtgN6.net
FPRPP.F <- function(r1,r2,n1,n2, prior=0.5){
p.value=fisher.test(matrix(c(r1,n1-r1,r2,n2-r2),2))$p.value
p1=r1/n1
p2=r2/n2
power=statmod::power.fisher.test(p1,p2,n1,n2,alpha=p.value)
FP=(1-prior)*p.value
TP=prior*power
c(p.value=p.value,'false positive rate' = FP/(TP+FP))
}
r1=10;n1=19;r2=11;n2=12
FPRPP.F(r1,r2,n1,n2,0.5)