17/10/25 19:11:29.66 nTDNW+TM.net
>>122
##
n=10
conf.level=0.95
PMF <- function(x) x^n # AUC==integrate(PMF,0,1) == 1/(n+1)
PDF <- function(x) (n+1)*x^n # PDF==PMF/AUC
Ex= (n+1)/(n+2) # integrate(x*PDF,0,1)==integrate((n+1)*x^(n+1),0,1)== (n+1)/(n+2)
CDF <- function(x) x^(n+1) # == integrate(PDF,0,x)
lwr = (1-conf.level)^(1/(n+1)) # CDF(lwr) == 1-conf.level
exp(log(1-conf.level)/(n+1))
curve(PDF)
curve(x*PDF(x))
integrate(function(x) x*PDF(x),0,1)$value
integrate(function(x) x*(n+1)*x^n,0,1)$value # integrate(function(x)(n+1)*x^(n+1),0,1)
(n+1)/(n+2)
# |((n+1)/(n+2))*x^(n+2)|[0,1]
201:卵の名無しさん
17/10/25 19:13:22.53 nTDNW+TM.net
>>193
## binomとの比較
binom::binom.confint(c(1,10,100),c(1,10,100),conf=0.95)
# n out of n
# mean=(n+1)/(n+2)
# lower=(n+1)√(1-conf.level) = 0.05^(1/(n+1))
non=function(n,conf.level=0.95){ # n hits out of n shoots
b=binom::binom.confint(n,n)[,c('method','mean','lower')]
n_1_root=function(n)(0.05)^(1/(n+1))
a=data.frame(method='root(n+1)',mean=(n+1)/(n+2),lower=(1-conf.level)^(1/(n+1)))
rbind(a,b)
}
> non(1)
method mean lower
1 root(n+1) 0.6666667 0.22360680
2 agresti-coull 1.0000000 0.16749949
3 asymptotic 1.0000000 1.00000000
4 bayes 0.7500000 0.22851981
5 cloglog 1.0000000 0.02500000
6 exact 1.0000000 0.02500000
7 logit 1.0000000 0.02500000
8 probit 1.0000000 0.02500000
9 profile 1.0000000 0.13811125
10 lrt 1.0000000 0.14650325
11 prop.test 1.0000000 0.05462076
12 wilson 1.0000000 0.20654931
202:卵の名無しさん
17/10/25 19:13:48.14 nTDNW+TM.net
> non(10)
method mean lower
1 root(n+1) 0.9166667 0.7615958
2 agresti-coull 1.0000000 0.6791127
3 asymptotic 1.0000000 1.0000000
4 bayes 0.9545455 0.8292269
5 cloglog 1.0000000 0.6915029
6 exact 1.0000000 0.6915029
7 logit 1.0000000 0.6915029
8 probit 1.0000000 0.6915029
9 profile 1.0000000 0.7303058
10 lrt 1.0000000 0.8252466
11 prop.test 1.0000000 0.6554628
12 wilson 1.0000000 0.7224672
> non(100)
method mean lower
1 root(n+1) 0.9901961 0.9707748
2 agresti-coull 1.0000000 0.9555879
3 asymptotic 1.0000000 1.0000000
4 bayes 0.9950495 0.9810231
5 cloglog 1.0000000 0.9637833
6 exact 1.0000000 0.9637833
7 logit 1.0000000 0.9637833
8 probit 1.0000000 0.9637833
9 profile 1.0000000 0.9670434
10 lrt 1.0000000 0.9809757
11 prop.test 1.0000000 0.9538987
12 wilson 1.0000000 0.9630065
203:卵の名無しさん
17/10/25 19:24:32.39 nTDNW+TM.net
>>195
n=10 の時のPDFと95%CI下限値0.7615958
URLリンク(i.imgur.com)
204:卵の名無しさん
17/10/25 19:24:57.43 nTDNW+TM.net
>>195
n=10 の時のPDFと95%CI下限値0.7615958
URLリンク(i.imgur.com)
205:卵の名無しさん
17/10/25 20:40:32.01 nTDNW+TM.net
#.n発r中の狙撃手が.N発狙撃するときの命中数を返す
Go13 <- function(.N, .n, r, k=10^3){ # k:シミュレーション回数
f <-function(S,N=.N,n=.n){
y=c(rep(1,S),rep(0,N-S))
sum(sample(y,n))
}
xx=r:.N
SS=NULL
for(i in 1:k){
x=sapply(xx,f)
SS=c(SS,which(x==r)-1+r)
}
print(summary(SS))
invisible(SS)
}
206:卵の名無しさん
17/10/26 12:27:43.35 ljbBi1cA.net
Golgo <- function(n=3,r=1,cl=0.95,K=10^6,Print=FALSE){
PMF <- function(x) choose(n,r)*x^r*(1-x)^(n-r)
AUC=integrate(PMF,0,1)$value
# library(hypergeo)
# f.auc <- function(x) choose(n,r)*x^(r+1)/(r+1)*hypergeo(r+1,r-n,r+2,x)
# auc=as.numeric(f.aunc(1)-f.auc(0)) ; 1/auc
PDF <- function(x)PMF(x)/AUC
Ex=integrate(function(x)x*PDF(x),0,1)$value
mode=optimize(PDF,c(0,1),maximum = TRUE)$maximum
CDF <- function(x) integrate(PDF,0,x)$value
CDFu0 <- function(x,u0=.025) CDF(x)-u0
lwr=uniroot(CDFu0,u0=(1-cl)/2,c(0,1))$root
upr=uniroot(CDFu0,u0=1-(1-cl)/2,c(0,1))$root
print(c(lower=lwr,mode=mode,mean=Ex,upper=upr))
xx=seq(0,1,by=0.0001)
pp=sapply(xx,PDF)
yy=sample(xx,K,replace=TRUE,prob=pp)
if(Print){hist(yy,freq=FALSE,main='',xlab='Probability',
col=sample(colors(),1),breaks=30)
lines(density(yy),col='gray',lty=3)
curve(PDF,add=TRUE)}
print(quantile(yy,probs=c(.025,.05,.5,.95,.975)))
myy=mean(yy)
dens <- density(yy)
mdyy=dens$x[which.max(dens$y)]
print(c(mode=mdyy,mean=myy))
invisible(yy)
}
207:卵の名無しさん
17/10/26 16:21:26.13 ljbBi1cA.net
# f回失敗した後にs回目の成功,K:シミュレーション回数
Dotsubo <- function(s=1,f=4,cl=0.95,K=10^6,Print=FALSE){
PMF <- function(p) choose(s+f-1,f)*(1-p)^f*p^s
AUC= (gamma(s+f)*gamma(s+1)) / (gamma(s)*gamma(f+s+2))
PDF <- function(p) (1-p)^f*p^s*gamma(f+s+2)/(gamma(f+1)*gamma(s+1))
curve(PDF)
Ex=integrate(function(x)x*PDF(x),0,1)$value
mode=optimize(PDF,c(0,1),maximum = TRUE)$maximum
CDF <- function(x) integrate(PDF,0,x)$value
CDFu0 <- function(x,u0=.025) CDF(x)-u0
lwr=uniroot(CDFu0,u0=(1-cl)/2,c(0,1))$root
upr=uniroot(CDFu0,u0=1-(1-cl)/2,c(0,1))$root
print(c(lower=lwr,mode=mode,mean=Ex,upper=upr))
xx=seq(0,1,by=0.0001)
pp=sapply(xx,PDF)
yy=sample(xx,K,replace=TRUE,prob=pp)
if(Print){hist(yy,freq=FALSE,main='',xlab='Probability',
col=sample(colors(),2),breaks=30)
lines(density(yy),col='gray',lty=3)
curve(PDF,add=TRUE)}
print(quantile(yy,probs=c(.025,.05,.5,.95,.975)))
myy=mean(yy)
dens <- density(yy)
mdyy=dens$x[which.max(dens$y)]
print(c(mode=mdyy,mean=myy))
invisible(yy)
}
208:卵の名無しさん
17/10/26 17:43:49.64 z4lFiDnk.net
ゴルゴ13は100発100中
ゴルゴ14は10発10中
ゴルゴ15は1発1中
とする。
各々10000発撃ったとき各ゴルゴの命中数の期待値はいくらか?
ドツボ13は100発0中
ドツボ14は10発0中
ドツボ15は1発0中
とする。
各々10000発撃ったときドツボの命中数の期待値はいくらか?
209:卵の名無しさん
17/10/27 15:54:12.90 dzxKDqmi.net
# n発r中の狙撃手がN発狙撃するときの命中数を返す
Golgo.sim <- function(N, n, r, k=10^3,Print=TRUE){ # k:シミュレーション回数
f <-function(S,.N=N,.n=n){ # 成績サンプル:命中S個、外れ(N-S)個
y=c(rep(1,S),rep(0,.N-S))
sum(sample(y,.n)) # その成績サンプルからn個数取り出したときの命中数
}
xx=r:N # r未満ではr個命中することはないので除外
SS=NULL # 容れ子
for(i in 1:k){
x=sapply(xx,f) # 命中数の配列
SS=append(SS,which(x==r)-1+r) # 命中数がr個のときの成績サンプルの命中数Sの配列をつくる
}
print(summary(SS))
print(quantile(SS,probs=c(0.025,0.05,0.95,0.975)))
print(c(mode=names(which.max(table(SS)))),quote=FALSE)
if(Print) {
hist(SS,xlim=c(0,N),freq=FALSE,col=sample(colors(),1),main='',xlab='Hits')
lines(density(SS))}
invisible(SS)
}
210:卵の名無しさん
17/10/27 15:55:45.80 pwjeiI6l.net
# n発r中の期待値
Golgo <- function(n=3,r=1,cl=0.95,K=10^6,Print=FALSE){
PMF <- function(x) choose(n,r)*x^r*(1-x)^(n-r)
AUC=integrate(PMF,0,1)$value
# library(hypergeo)
# f.auc <- function(x) choose(n,r)*x^(r+1)/(r+1)*hypergeo(r+1,r-n,r+2,x)
# auc=as.numeric(f.auc(1)-f.auc(0))
PDF <- function(x)PMF(x)/AUC
Ex=integrate(function(x)x*PDF(x),0,1)$value
mode=optimize(PDF,c(0,1),maximum = TRUE)$maximum
CDF <- function(x) integrate(PDF,0,x)$value
CDFu0 <- function(x,u0=.025) CDF(x)-u0
lwr=uniroot(CDFu0,u0=(1-cl)/2,c(0,1))$root
upr=uniroot(CDFu0,u0=1-(1-cl)/2,c(0,1))$root
print(c(lower=lwr,mode=mode,mean=Ex,upper=upr))
xx=seq(0,1,by=0.0001)
pp=sapply(xx,PDF)
yy=sample(xx,K,replace=TRUE,prob=pp)
if(Print){hist(yy,freq=FALSE,main='',xlab='Probability',
col=sample(colors(),1),breaks=30)
lines(density(yy),col='gray',lty=3)
curve(PDF,add=TRUE)}
print(quantile(yy,probs=c(.025,.05,.5,.95,.975)))
myy=mean(yy)
dens <- density(yy)
mdyy=dens$x[which.max(dens$y)]
print(c(mode=mdyy,mean=myy))
invisible(yy)
}
211:卵の名無しさん
17/10/29 17:45:50.97 LlbU36d2.net
data{// coin5.stan
int N;
int<lower=0,upper=1> Y[N];
}
parameters{
real<lower=0,upper=1> p;
}
model{
for(n in 1:N)
Y[n] ~ bernoulli(p);
}
data{// coin1.stan
int<lower=0,upper=1> Y;
}
parameters{
real<lower=0,upper=1> p;
}
model{
Y ~ bernoulli(p);
}
212:卵の名無しさん
17/10/29 22:16:10.95 LlbU36d2.net
# クラインの4元群Vが正6面体群S(P6)の正規部分群であることの確認
equal <- function(x,y) sum((x-y)^2)==0
tikan2 <- function(X,Y){ # replace Y first then X
n=length(X)
Z=rep(NA,n)
I=1:n
for(i in 1:n){
j <- which(I==X[i])
Z[i] <- Y[j]
}
return(Z)
}
f2 <- function(x,L=G){
re=NULL
for(i in 1:length(L)){
if(equal(L[[i]],x)) return(i) # i for names[i]
}
}
a =c(2,1,4,3)
b =c(3,4,1,2)
c =c(4,3,2,1)
e =c(1,2,3,4)
d =c(1,3,4,2)
t =c(1,2,4,3)
d2 =tikan2(d,d) # 1 4 2 3 # d3=tikan2(d2,d) == e
td =tikan2(t,d) # 1 3 2 4
td2=tikan2(td,d)
213:卵の名無しさん
17/10/29 22:17:23.29 LlbU36d2.net
V=list(e,a,b,c)
D3=list(e,d,d2,t,td,td2)
gr=expand.grid(V,D3)
.m=mapply(tikan2,gr[,1],gr[,2])
t.m=t(.m)
G=list()
for(i in 1:nrow(t.m)) G[[i]]=t.m[i,]
names(G)=paste0('t',1:length(G))
G # G: S(P6)
lG=length(G)
V
lV=length(V)
.M1=matrix(NA,nrow=lV,ncol=lG)
for (i in 1:lG){
for(j in 1:lV){
.M1[j,i]=f2(tikan2(G[[i]],V[[j]]))
}
}
.M2=matrix(NA,nrow=lV,ncol=lG)
for (i in 1:lG){
for(j in 1:lV){
.M2[j,i]=f2(tikan2(V[[j]],G[[i]]))
}
}
identical(apply(.M1,2,sort) ,apply(.M2,2,sort))
> identical(apply(.M1,2,sort) ,apply(.M2,2,sort))
[1] TRUE
214:卵の名無しさん
17/10/30 07:00:24.65 OjMfeTM6.net
## 写像を元とする集合の積を計算する
equal <- function(x,y) sum((x-y)^2)==0 & length(x)==length(y)
tikan2 <- function(X,Y){ # replace Y first then X
n=length(X)
Z=rep(NA,n)
I=1:n
for(i in 1:n){
j <- which(I==X[i])
Z[i] <- Y[j]
}
return(Z)
}
witch <- function(x,L=G){
re=NULL
for(i in 1:length(L)){
if(equal(L[[i]],x)) return(i) # i for names[i]
}
}
tikan3 <- function(X,Y,L=G){ # replace Y first then X and return index of L
Z=tikan2(X,Y)
witch(Z,L)
}
product <- function(A,B,L=G){ # product of sets, retur index of L
la=length(A)
lb=length(B
215:) gr=expand.grid(1:la,1:lb) f<-function(x,y) tikan3(A[[x]],B[[y]],L) re=mapply(f,gr[,1],gr[,2]) return(sort(unique(re))) }
216:卵の名無しさん
17/10/30 07:01:59.31 OjMfeTM6.net
GROUP <- function(n=6){
e=1:n
d1=c(n,1:(n-1)) # rotation
Mnd=matrix(NA,n,n)
Mnd[1,]=e
for(i0 in 1:(n-1)){
Mnd[i0+1,]=tikan(Mnd[i0,],d1)
}
rownames(Mnd)=c('e',paste0('d',1:(n-1)))
Mnd
t=n:1 ## mirror
tikan(t,d1) # D2n( e,d1,d2,d3,..,dn-1, t,td,td2,td3,...,tdn-1)
Mnt=matrix(NA,n,n)
Mnt[1,]=t
for(i2 in 1:(n-1)){
Mnt[i2+1,]=tikan(Mnt[i2,],d1)
}
Mnt
rownames(Mnt)=c('t',paste0('td',1:(n-1)))
Mn=rbind(Mnd,Mnt)
names=rownames(Mn)
Mn
.L=list(2*n) # list(Mn[1,],Mn[2,],...Mn[,n]) list(e,d1,d2,..,dn-1,t,td1,td2,..,tdn-1)
for(i3 in 1:(2*n)) .L[[i3]]=Mn[i3,]
names(.L)=names
G=.L
return(G)
}
217:卵の名無しさん
17/10/31 16:10:43.91 1m4iMpv8.net
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]))
}
218:卵の名無しさん
17/10/31 16:11:33.50 1m4iMpv8.net
# n発r中の期待値
Golgo2 <- function(n=3,r=1,cl=0.95,k=0.00001,Print=FALSE){
PMF <- function(x) choose(n,r)*x^r*(1-x)^(n-r)
xx=seq(0,1,by=k)
pmf=sapply(xx,PMF)
pdf=pmf/sum(pmf)
rspdf=rev(sort(pdf))
min.density=rspdf[min(which(cumsum(rspdf)>=cl))]
index=which(pdf>=min.density)
lower.idx=round(min(index))
upper.idx=round(max(index))
lower=xx[lower.idx]
upper=xx[upper.idx]
actual.CI=sum(pdf[index])
mpdf=sum(xx*pdf)
print(c(lower=lower,mean=mpdf,upper=upper,actual.CI=actual.CI),digits=5)
if(Print){plot(xx,pmf,lwd=1,xlab='probability of hit',
ylab=paste('probabilty of',r,'hits out of ',n,'shoots'))}
}
219:卵の名無しさん
17/10/31 20:38:26.69 A7oOP/mA.net
# 1年:進級失敗10人、うち1人放校
# 2年:進級失敗16人、放校なし
# 3年:進級失敗34人、うち放校9人
# 4年:進級失敗9人、うち放校2人
# 5年:進級失敗10人、うち放校1人
# 6年:卒業失敗26人、うち放校1人
# 一学年約120~130人前後。
#スレリンク(doctor板:1番)
flunk=c(10,16,34,9,10,26)
expel=c(1,0,9,2,1,1)
n=125
total.fl=sum(flunk)
total.ex=sum(expel)
re=matrix(NA,6,1)
for(i in 1:6){
re[i,]=poisson.test(c(flunk[i],n) ,c(total.fl,n*6))$p.value
}
re
for(i in 1:6){
re[i,]=poisson.test(c(expel[i],n) ,c(total.ex,n*6))$p.value
}
re
prop.test(flunk,rep(n,6))
pairwise.prop.test(flunk,rep(n,6),'bon')
prop.test(expel,rep(n,6))
d=cbind(expel,rep(n,6)-expel) ; d
fisher.test(d)
pairwise.prop.test(expel,rep(n,6),'none')
pairwise.prop.test(expel,rep(n,6),'bon')
220:卵の名無しさん
17/11/01 09:28:50.55 zVXhNw2e.net
## 確率分布から信頼区間を出す
HDCI2 <- function(PMF,cl=0.95,k=0.0001,Print=TRUE){
xx=seq(0,1,by=k)
xx=xx[-1]
pmf=sapply(xx,PMF)
pdf=pmf/sum(pmf)
rspdf=rev(sort(pdf))
min.density=rspdf[min(which(cumsum(rspdf)>=cl))]
index=which(pdf>=min.density)
lower.idx=round(min(index))
upper.idx=round(max(index))
lower=xx[lower.idx]
upper=xx[upper.idx]
actual.CI=sum(pdf[index])
mpdf=sum(xx*pdf)
print(c(lower=lower,mean=mpdf,upper=upper,actual.CI=actual.CI),digits=4)
if(Print) plot(xx,pmf,xlab='prior',ylab='posterior')
}
HDCI2(function(x)dnbinom(5-1,1,x))
221:卵の名無しさん
17/11/01 09:38:36.29 zVXhNw2e.net
## 確率分布から最尤値・期待値・信頼区間を出す
HDCI2 <- function(PMF,cl=0.95,k=0.0001,Print=FALSE){
pp=seq(0,1,by=k)
xx=pp[-1]
pmf=sapply(xx,PMF)
pdf=pmf/sum(pmf)
rspdf=rev(sort(pdf))
min.density=rspdf[min(which(cumsum(rspdf)>=cl))]
index=which(pdf>=min.density)
lower.idx=round(min(index))
upper.idx=round(max(index))
lower=xx[lower.idx]
upper=xx[upper.idx]
actual.CI=sum(pdf[index])
mpdf=sum(xx*pdf)
mode=xx[which.max(pmf)]
print(c(lower=lower,mode=mode,mean=mpdf,upper=upper,actual.CI=actual.CI),digits=4)
if(Print) plot(xx,pmf,xlab='prior',ylab='posterior')
}
222:卵の名無しさん
17/11/05 17:08:17.51 G4ZpDCNG.net
NHST Has 100% False Alarm Rate in Sequential
Testing
Under NHST, sequential testing of data generated from the null
hypothesis will eventually lead to a false alarm. With infinite
patience, there is 100% probability of falsely rejecting the null.
This is known as “sampling to reach a foregone conclusion” (e.g.,
Anscombe, 1954). To illustrate this phenomenon, a computer
simulation generated random values from a normal distribution
with mean zero and standard deviation one, assigning each sequential
value alternately to one or the other of two groups, and at each
step conducting a two-group t test assuming the current sample
sizes were fixed in advance. Each simulat
223:ed sequence began with N1 N2 3. If at any step the t test indicated p .05, the sequence was stopped and the total N ( N1 N2) was recorded. Otherwise, another random value was sampled from the zero-mean normal and included in the smaller group, and a new t test was conducted. For purposes of illustration, each sequence was limited to a maximum sample size of N 5,000. The simulation ran 1,000 sequences. NHST <- function(N=5000){ x=rnorm(3) y=rnorm(3) while(length(x) < N & t.test(x,y,var.equal=TRUE)$p.value >= 0.05){ x=append(x,rnorm(1)) y=append(y,rnorm(1)) } return(length(x)) } res <- replicate(1000,NHST())
224:卵の名無しさん
17/11/05 17:27:50.94 G4ZpDCNG.net
# NHST Has 100% False Alarm Rate in Sequential Testing (NHST : Null Hypothesis Significance Testing)
# Under NHST, sequential testing of data generated from the null
# hypothesis will eventually lead to a false alarm. With infinite
# patience, there is 100% probability of falsely rejecting the null.
# This is known as “sampling to reach a foregone conclusion” (e.g.,
# Anscombe, 1954). To illustrate this phenomenon, a computer
# a computer simulation generated random values from a normal distribution
# with mean zero and standard deviation one, assigning each sequential
# value alternately to one or the other of two groups, and at each
# step conducting a two-group t test assuming the current sample
# sizes were fixed in advance. Each simulated sequence began with
# N1=N2=3. If at any step the t test indicated p .05, the
# sequence was stopped and the total N (=N1+N2) was recorded.
NHST <- function(N){
x=rnorm(3)
y=rnorm(3)
while(length(x) < N & t.test(x,y,var.equal=TRUE)$p.value >= 0.05){
x=append(x,rnorm(1))
y=append(y,rnorm(1))
}
return(length(x))
}
FRR <- function(N,k=100){ #False Reject Rate
re <- replicate(k,NHST(N))
return(mean(re!=N)) # rate falsely rejected as significantly different
}
NN=c(10,20,40,80,160,320,640,1280,2560,5120,10240)
RjR=sapply(NN,FRR)
plot(NN,RjR,type='b')
225:卵の名無しさん
17/11/05 18:24:17.31 G4ZpDCNG.net
# Bayesian decision making, using the HDI and ROPE, does not
# suffer a 100% false alarm rate in sequential testing. Instead, the
# false alarm rate asymptotes at a much lower level, depending on
# the choice of ROPE. For illustration, again a computer simulation
# generated random values from a normal distribution with mean of
# zero and standard deviation of one, assigning each sequential value
# alternately to one or the other of two groups but at each step
# conducting a Bayesian analysis and checking whether the 95%
# HDI completely excluded or was contained within a ROPE from 0.15 to 0.15
# ROPE* region of practical equivalence
librry(BEST)
BA <- function(N){
mc=BEST::BESTmcmc(rnorm(N),rnorm(N),numSavedSteps=1000, burnInSteps=500)
re=summary(mc,ROPEm=c(-0.15,15))
return(re[['muDiff','%InROPE']])
}
NN=c(10,20,40,80,160,320,640,1280,2560,5120)
inrope=sapply(NN,BA)
plot(NN,inrope,type='b',pch=19,xlab='N',ylab='% in ROPE')
226:卵の名無しさん
17/11/05 20:18:58.81 G4ZpDCNG.net
Bayesian Estimation Supersedes the t-test (BEST) - online
URLリンク(www.sumsar.net)
トレースプロットがみられるのがうれしい。
227:卵の名無しさん
17/11/05 20:22:56.18 G4ZpDCNG.net
新版の予約受付中か。
URLリンク(www.amazon.co.jp)
ベイズ統計がどれくらい組み込まれたをみてから買うかな。
228:卵の名無しさん
17/11/05 20:28:14.21 G4ZpDCNG.net
# A君の彼女は女子大生、B君の彼女は女子高生。
# Y1女子大生n1=100人とY2女子高生n2=100人の胸囲を測定して
# 前者が平均82 , 標準偏差3
# 後者が平均81 , 標準偏差3
# 有意差はあるか?
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
2*pt(abs(T),n1-1+n2-1,lower.tail = FALSE)
}
T.test(100,100,82,81,3,3)
y1=82+scale(rnorm(100))*3
y2=81+scale(rnorm(100))*3
t.test(y1,y2,var.equal = TRUE)
library(BEST)
BESTout <- BESTmcmc(y1,y2)
plot(BESTout,ROPE=c(-2,2))
summary(BESTout,ROPEm=c(-2,2))
plotAll(BESTout,ROPEm=c(-2,2))
ProData <- makeData(mu1=82, sd1=3, mu2=81, sd2=3, nPerGrp=100,
pcntOut=10, sdOutMult=2.0, rnd.seed=NULL,showPlot=TRUE)
proMCMC <- BESTmcmc(proData$y1, proData$y2, numSavedSteps=2000)
N1plan <- N2plan <- 100
powerPro <- BESTpower(proMCMC, N1=N1plan, N2=N2plan,
ROPEm=c(-0.5,0.5), ROPEsd=c(-1,1), ROPEeff=c(-0.5,0.5),
maxHDIWm=5.0, maxHDIWsd=2.0, maxHDIWeff=1.0, nRep=5)
powerPro
229:卵の名無しさん
17/11/05 20:33:01.04 G4ZpDCNG.net
ProData <- makeData(mu1=82, sd1=3, mu2=81, sd2=3, nPerGrp=100,
pcntOut=10, sdOutMult=2.0, rnd.seed=NULL,showPlot=TRUE)
proMCMC <- BESTmcmc(proData$y1, proData$y2, numSavedSteps=2000)
N1plan <- N2plan <- 100
powerPro <- BESTpower(proMCMC, N1=N1plan, N2=N2plan,
ROPEm=c(-2,2), ROPEsd=c(-0,0), ROPEeff=c(-0,0),
maxHDIWm=5.0, maxHDIWsd=2.0, maxHDIWeff=1.0, nRep=1000)
powerPro
230:卵の名無しさん
17/11/05 20:51:07.47 G4ZpDCNG.net
Jikkan <- function(x){
re=summary(BESTout,ROPEm=c(-x,x))
re[['muDiff','%InROPE']]
}
xx=seq(0,2.5,by=0.01)
InRope=sapply(xx,Jikkan)
plot(xx,InRope,type='l',lwd=2, xlab='Difference',ylab='% Practically Equivalent')
abline(h=95, lty=3)
(Diff=uniroot(function(x,u0=95) Jikkan(x)-u0, c(1.5,2))$root)
plot(BESTout,ROPE=c(-Diff,Diff))
URLリンク(i.imgur.com)
231:卵の名無しさん
17/11/06 21:07:39.66 1uZMaFLW.net
##
library(BayesFactor)
tBF <- function(x,y){
t=t.test(x,y,var.equal = TRUE)$statistic
ttest.tstat(t,length(x),length(y),simple=TRUE)
}
NHST2 <- function(N){
x=rnorm(3)
y=rnorm(3)
while(length(x) < N & tBF(x,y) <= 1){
x=append(x,rnorm(1))
y=append(y,rnorm(1))
}
return(length(x))
}
FRR2 <- function(N,k=100){ #False Reject Rate
re <- replicate(k,NHST2(N))
return(mean(re!=N)) # rate falsely rejected as significantly different
}
> FRR2(10)
[1] 0.37
> FRR2(50)
[1] 0.46
> FRR2(100)
[1] 0.56
ベイズ因子をつかっても t 検定を組み込むと
# NHST Has 100% False Alarm Rate in Sequential Testing (NHST : Null Hypothesis Significance Testing)
という結論は変わらないな。
232:卵の名無しさん
17/11/06 21:51:42.88 1uZMaFLW.net
>>216
Sequential Samplingに修正
library(BEST)
BA <- function(x,y){
mc=BEST::BESTmcmc(x,y,num=5000)
re=summary(mc,ROPEm=c(-0.1,0.1))
return(re[['muDiff','%InROPE']])
}
x=rnorm(3)
y=rnorm(3)
inROPE=numeric()
N=100
for(i in 1:(N-3)){
x=append(x,rnorm(1))
y=append(y,rnorm(1))
inROPE[i]=BA(x,y)
}
plot(3:N,inROPE,type='b',pch=19,xlab='n',ylab='% in ROPE')
233:卵の名無しさん
17/11/07 00:05:34.55 2r/5Mjhu.net
>>201
n=1;r=1
f=function(p)choose(n,r)*p^r*(1-p)^(n-r)
AUC=integrate(f,0,1)$value
integrate(function(x) x*f(x)/AUC,0,1)$value
234:卵の名無しさん
17/11/07 00:49:20.56 2r/5Mjhu.net
library(BayesFactor)
tBF2 <- function(x,y,...){
tt=t.test(x,y,var.equal = TRUE)
bf=ttest.tstat(tt$statistic,length(x),length(y),simple=TRUE,...)
p=tt$p.value
return(c(bf=1/bf,p=p))
}
N=10
BF=numeric()
p.value=numeric()
k=1000
for(i in 1:k){
x=rnorm(N) ; y=rnorm(N)
bfp=tBF2(x,y)
BF[i]=bfp[1]
p.value[i]=bfp[2]
}
.m=cbind(BF,p.value)
.m[which(p.value<0.05),]
mean(p.value<0.05)
.m[which(BF<1),]
mean(BF<1)
235:卵の名無しさん
17/11/07 20:40:25.24 2r/5Mjhu.net
# crooked coin or dice
crooked <- function(n,r,H0=0.5,d=0.1,cl=0.95,Print=TRUE,...){
hdi=HDInterval::hdi(qbeta,cl,shape1=1+r,shape2=1+n-r)
if(Print){
ROPEl=H0*(1-d)
ROPEu=H0*(1+d)
curve(dbeta(x,1+r,1+n-r),lwd=1,xlab='p',ylab='density',...) # posterior
segments(ROPEl,0,ROPEu,lwd=5,col='navy')
segments(hdi[1],0,hdi[1],dbeta(hdi[1],1+r,1+n-r),col='blue')
segments(hdi[2],0,hdi[2],dbeta(hdi[2],1+r,1+n-r),col='blue')
legend('topleft',bty='n',legend=c('Range of Practival Equivalence',
'High Density Interval'),lwd=c(5,1),col=c('navy','blue'))
}
return(hdi)
}
> crooked(3000,490,H0=1/6,xlim=c(0.1,0.3))
lower upper
0.1503971 0.1768440
236:卵の名無しさん
17/11/07 21:15:32.89 2r/5Mjhu.net
# crooked coin or dice
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)
if(Print){
ROPEl=H0*(1-d)
ROPEu=H0*(1+d)
curve(dbeta(x,1+r,1+n-r),lwd=1,xlab='p',ylab='density',...) # posterior
segments(ROPEl,0,ROPEu,lwd=5,col='navy')
segments(hdi[1],0,hdi[1],dbeta(hdi[1],1+r,1+n-r),col='blue')
segments(hdi[2],0,hdi[2],dbeta(hdi[2],1+r,1+n-r),col='blue')
legend('topleft',bty='n',legend=c('Range of Practival Equivalence',
paste0(credMass/0.01,'% High Density Interval')),lwd=c(5,1),col=c('navy','blue'))
}
return(hdi)
}
crooked(3000,490,H0=1/6,xlim=c(0.12,0.20))
237:卵の名無しさん
17/11/08 19:29:03.43 vjtLqicz.net
data{// coin5d.stan
int N;
int<lower=0,upper=1> Y[N];
real b; // border
}
parameters{
real<lower=0,upper=1> p;
}
transformed parameters{
real ura = step(b-p); // ifelse(b-p>0,1,0)
}
model{
for(n in 1:N)
Y[n] ~ bernoulli(p);
}
238:卵の名無しさん
17/11/08 19:36:21.87 vjtLqicz.net
N=10
r=5
Y=c(rep(1,r),rep(0,N-r))
b=0.6
data <- list(Y=Y, N=N,b=b)
model.coin5d <- stan_model('coin5d.stan')
fit.coin5 <- sampling(model.coin5d, data=data, seed=123)
print(fit.coin5,digits=4)
Inference for Stan model: coin5d.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
p 0.4977 0.0035 0.1410 0.2286 0.3988 0.4966 0.5976 0.7744 1620 1.0036
ura 0.7552 0.0098 0.4300 0.0000 1.0000 1.0000 1.0000 1.0000 1933 1.0007
lp__ -8.8598 0.0234 0.7805 -11.0878 -9.0348 -8.5596 -8.3753 -8.3186 1116 1.0027
239:卵の名無しさん
17/11/08 20:57:34.40 vjtLqicz.net
# NHST Has 100% False Alarm Rate in Sequential Testing
# NHST : Null Hypothesis Significance Testing
# Under NHST, sequential testing of data generated from the null
# hypothesis will eventually lead to a false alarm. With infinite
# patience, there is 100% probability of falsely rejecting the null.
# This is known as “sampling to reach a foregone conclusion” (e.g.,
# Anscombe, 1954). To illustrate this phenomenon,
# a computer simulation generated random values from a normal distribution
# with mean zero and standard deviation one, assigning each sequential
# value alternately to one or the other of two groups, and at each
# step conducting a two-group t test assuming the current sample
# sizes were fixed in advance. Each simulated sequence began with
# N1=N2=3. If at any step the t test indicated p < .05, the
# sequence was stopped and the total N (=N1+N2) was recorded.
240:卵の名無しさん
17/11/08 20:57:42.57 vjtLqicz.net
FUN = wilcox.test
FUN = perm::permTS
FUN = lawstat::brunner.munzel.test
NHST <- function(N){
x=rnorm(3)
y=rnorm(3)
while(length(x) < N & FUN(x,y)$p.value >= 0.05){
x=append(x,rnorm(1))
y=append(y,rnorm(1))
}
return(length(x))
}
NHST(100)
FRR <- function(N,k=100){ #False Reject Rate
re <- replicate(k,NHST(N))
return(mean(re!=N)) # rate falsely rejected as significantly different
}
FRR(10)
FRR(50)
FRR(100)
FRR(500)
n=10 ; mean(replicate(k,FUN(rnorm(n),rnorm(n))$p.value) < 0.05)
n=50 ; mean(replicate(k,FUN(rnorm(n),rnorm(n))$p.value) < 0.05)
n=100 ; mean(replicate(k,FUN(rnorm(n),rnorm(n))$p.value) < 0.05)
n=500 ; mean(replicate(k,FUN(rnorm(n),rnorm(n))$p.value) < 0.05)
241:卵の名無しさん
17/11/11 09:46:29.46 EghxJpRr.net
source('URLリンク(github.com)')
source('URLリンク(github.com)')
a=1;b=1
N=21
z=21
pp=seq(0,1,by=0.001)
likelihood=pp^z*(1-pp)^(N-z)
prior=dbeta(pp,a,b)
posterior = BernGrid(pp,prior/sum(prior),c(rep(1,z),rep(0,N-z)), plotType="Bars" ,
showCentTend="Mean" , showHDI=TRUE , showpD=TRUE )
242:卵の名無しさん
17/11/11 10:09:38.03 EghxJpRr.net
# core concept
par(mfrow=c(2,2))
z=1 ; N=1
pp=seq(0,1,by=0.001)
prior=pmin(pp,1-pp) # a isosceles triangular distribution
plot(pp,prior,type='h',lwd=5,col='lightblue')
likelihood=pp^z*(1-pp)^(N-z)
plot(pp,likelihood,type='h',col='lightblue')
(pD=sum(prior*likelihood))
posterior=prior*likelihood/pD
plot(pp,posterior,type='h', col=' maroon')
243:卵の名無しさん
17/11/11 11:35:44.48 ARZ4w8db.net
require(rjags)
z=21
N=21
y=c(rep(1,z),rep(0,N-z))
data=list(N=N,y=y)
modelString = "
model {
for (i in 1:N) {
y[i] ~ dbern(theta)
}
theta ~ dbeta(1,1)
}
"
writeLines( modelString , con="TEMPmodel.txt" )
jagsModel = jags.model( file="TEMPmodel.txt" , data=data , # inits=initsList ,
n.chains=3 , n.adapt=500 )
update( jagsModel , n.iter=500 )
codaSamples = coda.samples( jagsModel , variable.names=c("theta") ,
n.iter=3334 )
summary(codaSamples)
plot(codaSamples)
244:卵の名無しさん
17/11/11 16:50:37.66 AQmX3Wf1.net
genMCMC <- # prior : theta ~ Beta(a,b)
function( data , numSavedSteps=50000 , saveName=NULL, a=1,b=1) {
require(rjags)
y = data$y
s = as.numeric(data$s) # converts character to consecutive integer levels
# Do some checking that data make sense:
if ( any( y!=0 & y!=1 ) ) { stop("All y values must be 0 or 1.") }
Ntotal = length(y)
Nsubj = length(unique(s))
# Specify the data in a list, for later shipment to JAGS:
dataList = list(
y = y ,
s = s ,
Ntotal = Ntotal ,
Nsubj = Nsubj
)
#-----------------------------------------------------------------------------
# THE MODEL.
modelString = paste0("
model {
for ( i in 1:Ntotal ) {
y[i] ~ dbern(
245: theta[s[i]] ) } for ( sIdx in 1:Nsubj ) { theta[sIdx] ~ dbeta(", a,',' , b," ) } } ") # close quote for modelString writeLines( modelString , con="TEMPmodel.txt" )
246:卵の名無しさん
17/11/11 16:50:57.13 AQmX3Wf1.net
#-----------------------------------------------------------------------------
# INTIALIZE THE CHAINS.
# Initial values of MCMC chains based on data:
# Option 1: Use single initial value for all chains:
# thetaInit = rep(0,Nsubj)
# for ( sIdx in 1:Nsubj ) { # for each subject
# includeRows = ( s == sIdx ) # identify rows of this subject
# yThisSubj = y[includeRows] # extract data of this subject
# thetaInit[sIdx] = sum(yThisSubj)/length(yThisSubj) # proportion
# }
# initsList = list( theta=thetaInit )
# Option 2: Use function that generates random values near MLE:
initsList = function() {
thetaInit = rep(0,Nsubj)
for ( sIdx in 1:Nsubj ) { # for each subject
includeRows = ( s == sIdx ) # identify rows of this subject
yThisSubj = y[includeRows] # extract data of this subject
resampledY = sample( yThisSubj , replace=TRUE ) # resample
thetaInit[sIdx] = sum(resampledY)/length(resampledY)
}
thetaInit = 0.001+0.998*thetaInit # keep away from 0,1
return( list( theta=thetaInit ) )
}
247:卵の名無しさん
17/11/11 16:51:08.92 AQmX3Wf1.net
#-----------------------------------------------------------------------------
# RUN THE CHAINS
parameters = c( "theta") # The parameters to be monitored
adaptSteps = 500 # Number of steps to adapt the samplers
burnInSteps = 500 # Number of steps to burn-in the chains
nChains = 4 # nChains should be 2 or more for diagnostics
thinSteps = 1
nIter = ceiling( ( numSavedSteps * thinSteps ) / nChains )
# Create, initialize, and adapt the model:
jagsModel = jags.model( "TEMPmodel.txt" , data=dataList , inits=initsList ,
n.chains=nChains , n.adapt=adaptSteps )
# Burn-in:
cat( "Burning in the MCMC chain...\n" )
update( jagsModel , n.iter=burnInSteps )
# The saved MCMC chain:
cat( "Sampling final MCMC chain...\n" )
codaSamples = coda.samples( jagsModel , variable.names=parameters ,
n.iter=nIter , thin=thinSteps )
# resulting codaSamples object has these indices:
# codaSamples[[ chainIdx ]][ stepIdx , paramIdx ]
if ( !is.null(saveName) ) {
save( codaSamples , file=paste(saveName,"Mcmc.Rdata",sep="") )
}
return( codaSamples )
}
248:卵の名無しさん
17/11/11 19:06:46.40 AQmX3Wf1.net
library("rstan")
library("rstudioapi")
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
G1314model='
data{// Golgo13.14.stan
int N13; // 100
int N14; //10
int<lower=0,upper=1> G13[N13];
int<lower=1,upper=2> G14[N14];
}
parameters{
real<lower=0,upper=1> p13;
real<lower=0,upper=1> p14;
}
model {
for(i in 1:N13) G13[i] ~ bernoulli(p13);
for(j in 1:N14) G14[j] ~ bernoulli(p14);
p13 ~ beta(1,1);
p14 ~ beta(1,1);
}
generated quantities{
real d = p13 - p14;
}
'
writeLines( G1314model , con='G1314.stan' )
249:卵の名無しさん
17/11/11 19:15:36.64 AQmX3Wf1.net
N13=100
N14=10
G13=rep(1,N13)
G14=rep(1,N14)
data=list(N13=N13,N14=N14,G13=G13,G14=G14)
G1314.model=stan_model('G1314.stan')
saveRDS(G1314.model,file='G1314.model.rds')
# G1314.model=readRDS('G1314.model.rds')
fitG1314 <- sampling(G1314.model, data=data, seed=1234)
fitG1314
stan_dens(fitG1314,separate_chains = TRUE)
i.imgur.com/aL2PSDM.png
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
p13 0.99 0.00 0.01 0.97 0.99 0.99 1.00 1.00 2500 1
p14 0.92 0.00 0.08 0.72 0.88 0.94 0.97 1.00 2834 1
d 0.07 0.00 0.08 -0.01 0.02 0.05 0.11 0.27 2837 1
lp__ -10.24 0.04 1.16 -13.44 -10.70 -9.88 -9.42 -9.09 979 1
> ms=rstan::extract(fitG1314)
> mean(ms$d<0)
[1] 0.09425
250:卵の名無しさん
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
1次方程式もできないド底辺特殊シリツ医大卒の記録
URLリンク(imagizer.imageshack.com)
何度読んでも馬鹿すぎる。
男女別の割合と全体での割合から男女比が計算できるとも思わないとは。
なんでこんなのが大学に入れるわけよ?
裏口入学以外に説明がつく?
中学生でも解ける一次方程式の問題だろ。
それすらできない馬鹿が自信を持って発言。
>患者の男女比が必要なのもわからないのか?
だとさ。
URLリンク(imagizer.imageshack.com)
0.2575549
と答を書いてやったら
>単位も書かずに答えだとか…
ド底辺シリツ医大では確率に単位があるらしいぞwww
何でこんな馬鹿が大学に入れるわけ?
裏口入学以外に説明がつく?
URLリンク(imagizer.imageshack.com)
298:卵の名無しさん
17/12/03 08:22:35.16 Egt6Q5KK.net
1次方程式もできないド底辺特殊シリツ医大卒の記録
URLリンク(imagizer.imageshack.com)
何度読んでも馬鹿すぎる。
男女別の割合と全体での割合から男女比が計算できるとも思わないとは。
なんでこんなのが大学に入れるわけよ?
裏口入学以外に説明がつく?
中学生でも解ける一次方程式の問題だろ。
シリツ医大には二次方程式が解けないやつがいると言ってた えなりかずき もビックリだろね。
それすらできない馬鹿が自信を持って発言。
>患者の男女比が必要なのもわからないのか?
だとさ。
URLリンク(imagizer.imageshack.com)
求める確率を
0.2575549
と答を書いてやったら
>単位も書かずに答えだとか…
ド底辺シリツ医大では確率に単位があるらしいぞwww
何でこんな馬鹿が大学に入れるわけ?
裏口入学以外に説明がつく?
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
5%幅の違いは同等扱いにすると
> 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
と0を挟む。
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
5%幅の違いは同等扱いにすると
> 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
と0を挟む。
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
5%幅の違いは同等扱いにすると
> 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
と0を挟む。
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
平均・標準偏差が以下の4群の多重比較
> 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
p値で考えないから補正も不要のベイズの方が直感的だな。
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が0を跨いでしまう。
URLリンク(i.imgur.com)
316:卵の名無しさん
17/12/16 16:52:42.18 LibxCzfo.net
グループごとに固有の標準偏差を持つとすると
URLリンク(i.imgur.com)
で HDIが0を跨がない。
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)
べつに分布を90度回転させて表示させなくてもいいのだが。
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
# 薬剤yは3人目で治癒、薬剤gは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人中z人開脚、どっちが開脚が容易か?
という問題にした方が興味をひくなw
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 人中 z 人が開脚、どっちが開脚が容易か?
r=3
z=3
N=9
とサンプルでの比率が同じとき母集団の推定平均値に差があるだろうか?
stanの出力をグラフにしてみた。
平均値で4%ほどの差が推定された。
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
あるド底辺シリツ医大で入学者の裏口入学者と学力考査合格入学者の比率は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
で帰無仮説は棄却できないと結論した。
URLリンク(i.imgur.com)
一方、本番と十八番が好きな太郎は一人ずつ調べて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
URLリンク(i.imgur.com)
で帰無仮説は棄却される。
331:卵の名無しさん
17/12/20 21:16:42.01 e+3oE/TR.net
コインが続けて5回裏がでたときにこのコインはイカサマコインといえるか?
5%のばらつきはイカサマとみなさないとする。
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の(合格者数+1)乗根
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)回コインをなげてn(=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
サイコロの1の目が何回続けてでたらイカサマサイコロかを考えて暇つぶし。
このスレの趣旨wに合わせてこんな問題にしてみた。
ド底辺学力学生がド底辺特殊シリツ医大を受験したとする。
試験は五者択一で三教科150問、合格ラインは6割とされる。
中学数学すらできないド底辺学力ゆえ正解できるのは偶然に頼るしか�
347:ネく、 正答率は概ね1/5で、その日の運で1/4から1/6と推定されている。 これを、正答確率は最頻値1/5で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
5試合連続で勝敗予想的中なら頻度主義では予知能力あるとされる。p=0.03125 < 0.05
URLリンク(to-kei.net)
ベイズでやってみるなら
的中率が1/2である確率は一応分布に従う(事前分布)として
5試合連続的中した後の的中率事後分布がどうなるかを考える。
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')