DBDA with pymc3

from IPython.core.pylabtools import figsize
import numpy as np
from matplotlib import pyplot as plt
from os import makedirs
from urllib.request import urlretrieve
import pymc3 as pm
import theano.tensor as tt
from pymc3 import model_to_graphviz
from pymc3 import forestplot


## analysis 1 P208
urlretrieve("http://peach.l.chiba-u.ac.jp/course_folder/z15N50.csv", "data/z15N50.csv")
data01 = np.loadtxt("data/z15N50.csv", skiprows=1)
print(data01)

with pm.Model() as model:
    theta = pm.Beta('theta', alpha = 1, beta = 1)
    obs = pm.Bernoulli('obs', theta, observed = data01)

model_to_graphviz(model)

with model:
    trace = pm.sample(10000)


# checking results
figsize(12.5, 4)
pm.plots.traceplot(trace, var_names=["theta"])
pm.plots.plot_posterior(trace, var_names=["theta"])
pm.plots.autocorrplot(trace, var_names=["theta"])


# Analysis 2 p208
resp = [1,0,1,1,1,1,1,0,0,0,1,0,0,1,0]
subj = [0,0,0,0,0,0,0,0,1,1,1,1,1,1,1]

from pymc3 import model_to_graphviz
N_subj=2
with pm.Model() as model2:

    theta = pm.Beta('theta', alpha = 2, beta = 2, shape=N_subj)

    y = pm.Bernoulli('y', theta,  observed=resp)
    
    diff_theta = pm.Deterministic('delta_theta', theta[0] - theta[1])
model_to_graphviz(model2)

# mcmc sampling
with model2:
    trace = pm.sample(10000)

# viz & checking results
pm.plots.traceplot(trace, var_names=["theta"])
pm.plots.plot_posterior(trace, var_names=["theta"])
pm.plots.plot_posterior(trace, var_names=["delta_theta"])
pm.plots.autocorrplot(trace, var_names=["theta"])


### analysis 3 p240
N_subj = len(data.s.unique())
resp = data.y
subj = data.s - 1
with pm.Model() as model3:
    omega = pm.Beta('omega',1,1)
    kappaM2 = pm.Gamma('kappaM2', 0.01, 0.01)
    kappa = pm.Deterministic('kappa', kappaM2 + 2)
    theta = pm.Beta('theta', 
                    alpha = omega*(kappa-2)+1,
                    beta = (1-omega)*(kappa-2)+1, 
                    shape = N_subj)



    y = pm.Bernoulli('y', theta[subj],  observed=resp)
    
model_to_graphviz(model3)
### analysis 4 p253
N_player = len(data.Player.unique())
N_pos = 9
Ns = data.AtBats
Zs = data.Hits
subj = data.PlayerNumber - 1
pos = data.PriPosNumber - 1

with pm.Model() as model4:
    omegaO = pm.Beta('omegaO',1,1)
    kappaM2O = pm.Gamma('kappaM2', 0.01, 0.01)
    kappaO = pm.Deterministic('kappaO', kappaM2O + 2)
    
    omegaPOS = pm.Beta('omegaPOS', 
                       alpha = omegaO*(kappaO-2)+1, 
                       beta = (1-omegaO)*(kappaO-2)+1, 
                       shape = N_pos)
    kappaM2POS = pm.Gamma('kappaM2POS', 0.01, 0.01,
                           shape = N_pos)
    kappaPOS = pm.Deterministic('kappaPOS', kappaM2POS + 2)

   
    theta = pm.Beta('theta', 
                    alpha = omegaPOS[pos]*(kappaPOS[pos]-2)+1, 
                   beta = (1-omegaPOS[pos])*(kappaPOS[pos]-2)+1, 
                    shape = N_player)

    y = pm.Binomial('y', n = Ns, p = theta,  observed = Zs)
    
model_to_graphviz(model4)

院:認識情報解析

library(rjags)
source("http://peach.l.chiba-u.ac.jp/course_folder/HDI_revised.txt")

dat<-read.csv("http://www.matsuka.info/data_folder/HtWtData110.csv")
library(plot3D)
w = seq(80,360,length.out=100)
h = seq(50, 75, length.out=100)
M <- mesh(w,h)
P.male = 1/(1+exp(-1*(0.018*M$x+0.7*M$y-50)))

scatter3D(dat$weight, dat$height, dat$mal, pch = 19, cex = 2,
          theta = 30, phi = 45, ticktype = "detailed", zlim=c(-0.1,1),ylim=c(50,78),xlim=c(80,360),
          xlab = "weight", ylab = "height", zlab = "P(male)",
          surf = list(x = M$x, y = M$y, z = P.male,facets = NA))


y = dat$male; x = dat$weight; Ntotal = length(y)
dataList = list(y = y, x = x, Ntotal = Ntotal)

model.txt = "
data {
  xm <- mean(x)
  xsd <- sd(x)
  for (i in 1:Ntotal){
    zx[i] = (x[i] - xm)/xsd
  }
}
model {
  for ( i_data in 1:Ntotal ) {
    y[ i_data ] ~ dbern(ilogit( zbeta0 + zbeta * zx[i_data]))
  }
  zbeta0 ~ dnorm(0, 1/2^2)
  zbeta ~ dnorm(0, 1/2^2)

  beta <- zbeta / xsd
  beta0 <- zbeta0 - zbeta * xm/xsd
}"
writeLines(model.txt, "model1.txt")
parameters = c( "beta0" ,  "beta")
jagsModel = jags.model( "model1.txt", data=dataList, n.chains=3, n.adapt=500 )
update( jagsModel , n.iter=1000)
codaSamples = coda.samples( jagsModel , variable.names=parameters, n.iter=10000, thin=5)
mcmcMat<-as.matrix(codaSamples)

plot(dat$weight,dat$male,xlim=c(90,280),yaxt="n",ylab="Male / Female",
     xlab="Weight", cex=2.5)
axis(2,at = 0:1,labels=c("Femal","Male"))
n2plot=100
idx=sample(1:nrow(mcmcMat),n2plot)
temp.x = seq(90,280,length.out = 100)
for (i_sample in 1:n2plot) {
  temp.y = 1/(1+exp(-1*(mcmcMat[idx[i_sample],2] + mcmcMat[idx[i_sample],1]*temp.x)))
  lines(temp.x, temp.y, col='orange', lwd=2)
}


x = cbind(dat$weight,dat$height);Nx = ncol(x)
dataList = list(y = y, x = x, Ntotal = Ntotal, Nx = Nx)

model.txt = "
data {
  for (j in 1:Nx){
    xm[j] <- mean(x[,j])
    xsd[j] <- sd(x[,j])
    for (i in 1:Ntotal){
      zx[i,j] = (x[i,j] - xm[j])/xsd[j]
    }
  }
}
model {
  for ( i_data in 1:Ntotal ) {
    y[ i_data ] ~ dbern(ilogit( zbeta0 + sum(zbeta[1:Nx] * zx[i_data, 1:Nx ])))
  }
  zbeta0 ~ dnorm(0, 1/2^2)
  for (j in 1:Nx){
    zbeta[j] ~ dnorm(0, 1/2^2)
  }
  beta[1:Nx] <- zbeta[1:Nx] / xsd[1:Nx]
  beta0 <- zbeta0 -sum(zbeta[1:Nx] * xm[1:Nx]/xsd[1:Nx])
}"
writeLines(model.txt, "model.txt")

parameters = c( "beta0" ,  "beta")
jagsModel = jags.model( "model.txt", data=dataList, n.chains=3, n.adapt=500 )
update( jagsModel , n.iter=1000)
codaSamples = coda.samples( jagsModel , variable.names=parameters, n.iter=10000, thin=5)
mcmcMat<-as.matrix(codaSamples)
par(mfrow=c(1,3))
HDI.plot(mcmcMat[,3],xlabel='intercept')
HDI.plot(mcmcMat[,1],xlabel='weight')
HDI.plot(mcmcMat[,2],xlabel='height')

par(mfrow=c(1,1))
plot(dat$weight,dat$height,xlab="Weight", ylab="Height", type="n")
n2plot=100
idx=sample(1:nrow(mcmcMat),n2plot)
for (i_sample in 1:n2plot) {
  abline(a=-1*mcmcMat[idx[i_sample],3]/mcmcMat[idx[i_sample],2],
         b=-1*mcmcMat[idx[i_sample],1]/mcmcMat[idx[i_sample],2],col="orange")

}
points(dat$weight,dat$height,pch=paste(dat$male), cex=1.5)

# un-even data
x = rnorm(300)
pr = 1/(1+exp(2*x))
y = pr < runif(300)
plot(x,y)

remove.id = sample(which(y == 0),120)

Ntotal = length(y[-remove.id])
dataList = list(y = y[-remove.id], x = x[-remove.id], Ntotal = Ntotal)

model.txt = "
data {
xm <- mean(x)
xsd <- sd(x)
for (i in 1:Ntotal){
zx[i] = (x[i] - xm)/xsd
}
}
model {
for ( i_data in 1:Ntotal ) {
y[ i_data ] ~ dbern(ilogit( zbeta0 + zbeta * zx[i_data]))
}
zbeta0 ~ dnorm(0, 1/2^2)
zbeta ~ dnorm(0, 1/2^2)

beta <- zbeta / xsd
beta0 <- zbeta0 - zbeta * xm/xsd
}"
writeLines(model.txt, "model1.txt")
parameters = c( "beta0" ,  "beta")
jagsModel = jags.model( "model1.txt", data=dataList, n.chains=3, n.adapt=500 )
update( jagsModel , n.iter=1000)
codaSamples = coda.samples( jagsModel , variable.names=parameters, n.iter=10000, thin=5)
mcmcMat<-as.matrix(codaSamples)

plot(x[-remove.id],y[-remove.id],xlim=c(-3,3))
n2plot=100
idx=sample(1:nrow(mcmcMat),n2plot)
temp.x = seq(-3,3,length.out = 100)
for (i_sample in 1:n2plot) {
  temp.y = 1/(1+exp(-1*(mcmcMat[idx[i_sample],2] + mcmcMat[idx[i_sample],1]*temp.x)))
  lines(temp.x, temp.y, col='orange', lwd=2)
}


x1 = rnorm(150)
x2 = x1*0.9+rnorm(150,0,0.5)
pr = 1/(1+exp(x1+x2))
y = pr < runif(150)
Ntotal = length(y)
dataList = list(y = y, x = cbind(x1,x2), Ntotal = Ntotal, Nx = 2)


model.txt = "
data {
for (j in 1:Nx){
xm[j] <- mean(x[,j])
xsd[j] <- sd(x[,j])
for (i in 1:Ntotal){
zx[i,j] = (x[i,j] - xm[j])/xsd[j]
}
}
}
model {
for ( i_data in 1:Ntotal ) {
y[ i_data ] ~ dbern(ilogit( zbeta0 + sum(zbeta[1:Nx] * zx[i_data, 1:Nx ])))
}
zbeta0 ~ dnorm(0, 1/2^2)
for (j in 1:Nx){
zbeta[j] ~ dnorm(0, 1/2^2)
}
beta[1:Nx] <- zbeta[1:Nx] / xsd[1:Nx]
beta0 <- zbeta0 -sum(zbeta[1:Nx] * xm[1:Nx]/xsd[1:Nx])
}"
writeLines(model.txt, "model.txt")

parameters = c( "beta0" ,  "beta")
jagsModel = jags.model( "model.txt", data=dataList, n.chains=3, n.adapt=500 )
update( jagsModel , n.iter=1000)
codaSamples = coda.samples( jagsModel , variable.names=parameters, n.iter=10000, thin=5)
mcmcMat<-as.matrix(codaSamples)
plot(x1,x2,xlab="x1", ylab="x2", type="n")
n2plot=100
idx=sample(1:nrow(mcmcMat),n2plot)
for (i_sample in 1:n2plot) {
  abline(a=-1*mcmcMat[idx[i_sample],3]/mcmcMat[idx[i_sample],2],
         b=-1*mcmcMat[idx[i_sample],1]/mcmcMat[idx[i_sample],2],col="orange")

}
points(x1,x2,pch=paste(y), cex=1.5)

# guessing
y = dat$male; x = dat$weight; Ntotal = length(y)
dataList = list(y = y, x = x, Ntotal = Ntotal)

model.txt = "
data {
xm <- mean(x)
xsd <- sd(x)
for (i in 1:Ntotal){
zx[i] = (x[i] - xm)/xsd
}
}
model {
for ( i_data in 1:Ntotal ) {
y[ i_data ] ~ dbern(mu[i_data])
mu[i_data] <- (guess*0.5 + (1-guess)*ilogit( zbeta0 + zbeta * zx[i_data]))
}
zbeta0 ~ dnorm(0, 1/2^2)
zbeta ~ dnorm(0, 1/2^2)
guess ~ dbeta(1,9)
beta <- zbeta / xsd
beta0 <- zbeta0 - zbeta * xm/xsd
}"
writeLines(model.txt, "model1.txt")
parameters = c( "beta0" ,  "beta", "guess")
jagsModel = jags.model( "model1.txt", data=dataList, n.chains=3, n.adapt=500 )
update( jagsModel , n.iter=1000)
codaSamples = coda.samples( jagsModel , variable.names=parameters, n.iter=10000, thin=5)
mcmcMat<-as.matrix(codaSamples)

plot(x[-remove.id],y[-remove.id],xlim=c(-3,3))
n2plot=100
idx=sample(1:nrow(mcmcMat),n2plot)
temp.x = seq(-3,3,length.out = 100)
for (i_sample in 1:n2plot) {
  temp.y = 1/(1+exp(-1*(mcmcMat[idx[i_sample],2] + mcmcMat[idx[i_sample],1]*temp.x)))
  lines(temp.x, temp.y, col='orange', lwd=2)
}
par(mfrow=c(1,3))
HDI.plot(mcmcMat[,2],xlabel='intercept')
HDI.plot(mcmcMat[,1],xlabel='weight')
HDI.plot(mcmcMat[,3],xlabel='guessing')

par(mfrow=c(1,1))
plot(dat$weight,dat$male,xlim=c(90,280),yaxt="n",ylab="Male / Female",
     xlab="Weight", cex=2.5)
axis(2,at = 0:1,labels=c("Femal","Male"))
n2plot=100
idx=sample(1:nrow(mcmcMat),n2plot)
temp.x = seq(90,280,length.out = 100)
for (i_sample in 1:n2plot) {
  temp.y = mcmcMat[idx[i_sample],3]/2+(1-mcmcMat[idx[i_sample],3])*1/(1+exp(-1*(mcmcMat[idx[i_sample],2] + mcmcMat[idx[i_sample],1]*temp.x)))
  lines(temp.x, temp.y, col='orange', lwd=2)
}

# nomial predictors
model.txt = "
model {
  for ( i.data in 1:Ntotal ) {
    y[ i.data ] ~ dbin(mu[i.data],N[i.data])
    mu[i.data] ~ dbeta(omega[x[i.data]]*(kappa-2)+1,(1-omega[x[i.data]])*(kappa-2)+1)
  }
  for (i.pos in 1:Npos){
    omega[i.pos] <- ilogit(a0+a[i.pos])
    a[i.pos] ~ dnorm(0.0, 1/aSigma^2)
  }
  a0 ~  dnorm(0,1/2^2)
  aSigma ~ dgamma(1.64, 0.32)
  kappa <- kappaMinusTwo +2
  kappaMinusTwo ~ dgamma(0.01,0.01)
  for (i.pos in 1:Npos){
    m[i.pos] <- a0+a[i.pos]
  }
  b0 <- mean(m[1:Npos])
  for (i.pos in 1:Npos){
    b[i.pos] <- m[i.pos] - b0
  }
}"
writeLines(model.txt, "model.txt")

dat<-read.csv("http://www.matsuka.info/data_folder/BattingAverage.csv")
y = dat$Hits
N = dat$AtBats
x = dat$PriPosNumber
Ntotal = length(y)
Npos = length(unique(x))
dataList = list(y = y, x = x, N = N, Ntotal = Ntotal, Npos = Npos)
parameters = c( "b0" ,  "b", "omega")
jagsModel = jags.model( "model.txt", data=dataList, n.chains=3, n.adapt=500 )
update( jagsModel , n.iter=1000)
codaSamples = coda.samples( jagsModel , variable.names=parameters, n.iter=10000, thin=5)
mcmcMat<-as.matrix(codaSamples)

par(mfrow=c(3,3))
for (i.pos in 1:9){
  HDI.plot(mcmcMat[,i.pos+10])
}

par(mfrow=c(2,2))
HDI.plot(mcmcMat[,1]-mcmcMat[,2])
HDI.plot(mcmcMat[,2]-mcmcMat[,3])
HDI.plot(mcmcMat[,11]-mcmcMat[,12])
HDI.plot(mcmcMat[,12]-mcmcMat[,13])

# softmax regression
x1 = runif(500, min=-2, max = 2)
x2 = runif(500, min=-2, max = 2)
b0 = c(0,-3,-4,-5)
b1 = c(0,-5,-1,10)
b2 = c(0,-5,10,-1)
l1 = b0[1]+b1[1]*x1+b2[1]*x2
l2 = b0[2]+b1[2]*x1+b2[2]*x2
l3 = b0[3]+b1[3]*x1+b2[3]*x2
l4 = b0[4]+b1[4]*x1+b2[4]*x2
p1 = exp(l1)/sum(exp(l1)+exp(l2)+exp(l3)+exp(l4))
p2 = exp(l2)/sum(exp(l1)+exp(l2)+exp(l3)+exp(l4))
p3 = exp(l3)/sum(exp(l1)+exp(l2)+exp(l3)+exp(l4))
p4 = exp(l4)/sum(exp(l1)+exp(l2)+exp(l3)+exp(l4))
ps = cbind(p1,p2,p3,p4)
y = apply(ps,1,which.max)
plot(x1,x2,pch=y,col=y)

b0 = c(0,-4,-1,-1)
b1 = c(0,-5,1,3)
b2 = c(0,0,-5,3)
l1 = b0[1]+b1[1]*x1+b2[1]*x2
l2 = b0[2]+b1[2]*x1+b2[2]*x2
l3 = b0[3]+b1[3]*x1+b2[3]*x2
l4 = b0[4]+b1[4]*x1+b2[4]*x2
p1 = exp(l1)/sum(exp(l1)+exp(l2)+exp(l3)+exp(l4))
p2 = exp(l2)/sum(exp(l1)+exp(l2)+exp(l3)+exp(l4))
p3 = exp(l3)/sum(exp(l1)+exp(l2)+exp(l3)+exp(l4))
p4 = exp(l4)/sum(exp(l1)+exp(l2)+exp(l3)+exp(l4))
ps = cbind(p1,p2,p3,p4)
y = apply(ps,1,which.max)
plot(x1,x2,pch=y,col=y)

p1 = exp(l1)/sum(exp(l1)+exp(l2)+exp(l3)+exp(l4))
p2 = exp(l2)/sum(exp(l1)+exp(l2)+exp(l3)+exp(l4))
p3 = exp(l3)/sum(exp(l1)+exp(l2)+exp(l3)+exp(l4))
p4 = exp(l4)/sum(exp(l1)+exp(l2)+exp(l3)+exp(l4))
ps = cbind(p1,p2,p3,p4)
p12 = pmax(p1,p2)
p34 = pmax(p3,p4)
y12vs34 = apply(cbind(p1,p2),1,which.max)
plot(x1,x2,pch=y12vs34,col=y12vs34)
y1vs2 = apply(cbind(p1,p3),1,which.max)
points(x1,x2,pch=y1vs2+2,col=y1vs2+2)
y3vs4 = apply(cbind(p1,p4),1,which.max)
points(x1,x2,pch=y3vs4+6,col=y3vs4+6)



model.txt = "
data {
  for ( j in 1:Nx ) {
    xm[j]  <- mean(x[,j])
    xsd[j] <-   sd(x[,j])
    for ( i in 1:Ntotal ) {
      zx[i,j] <- ( x[i,j] - xm[j] ) / xsd[j]
    }
  }
}
model {
  for ( i in 1:Ntotal ) {
    y[i] ~ dcat(mu[1:Nout,i])
    mu[1:Nout,i] <- explambda[1:Nout,i]/sum(explambda[1:Nout,i])
    for (k in 1:Nout){
      explambda[k,i]=exp(zbeta0[k] + sum(zbeta[k,1:Nx] * zx[i, 1:Nx ]))
    }
  }
  zbeta0[1] = 0
  for (j in 1:Nx){
    zbeta[1,j] <- 0
  }
  for (k in 2:Nout){
    zbeta0[k] ~ dnorm(0, 1/2^2)
    for (j in 1:Nx){
      zbeta[k,j]~dnorm(0, 1/2^2)
    }
  }
  for ( k in 1:Nout ) {
    beta[k,1:Nx] <- zbeta[k,1:Nx] / xsd[1:Nx]
    beta0[k] <- zbeta0[k] - sum( zbeta[k,1:Nx] * xm[1:Nx] / xsd[1:Nx] )
  }
}"
writeLines(model.txt, "model.txt")

dat<-read.csv( "http://www.matsuka.info/data_folder/CondLogistRegData1.csv" )
y = dat$Y
x = cbind(dat[,1],dat[,2])
Ntotal = length(y)
Nout = length(unique(y))
dataList = list(y = y, x = x, Nx = 2, Ntotal = Ntotal, Nout = Nout)
parameters = c( "beta0" ,  "beta")
jagsModel = jags.model( "model.txt", data=dataList, n.chains=3, n.adapt=500 )
update( jagsModel , n.iter=1000)
codaSamples = coda.samples( jagsModel , variable.names=parameters, n.iter=10000, thin=5)
mcmcMat<-as.matrix(codaSamples)
par(mfrow=c(1,3))
HDI.plot(mcmcMat[,7+0],xlab='intercept')
HDI.plot(mcmcMat[,1+0],xlab='b1')
HDI.plot(mcmcMat[,4+0],xlab='b2')

model = "
data {
  for ( j in 1:Nx ) {
    xm[j]  <- mean(x[,j])
    xsd[j] <-   sd(x[,j])
    for ( i in 1:Ntotal ) {
      zx[i,j] <- ( x[i,j] - xm[j] ) / xsd[j]
    }
  }
}
# Specify the model for standardized data:
model {
  for ( i in 1:Ntotal ) {
    y[i] ~ dcat( mu[1:Nout,i] )
    mu[1,i] <- phi[1,i]
    mu[2,i] <- phi[2,i] * (1-phi[1,i])
    mu[3,i] <- phi[3,i] * (1-phi[2,i]) * (1-phi[1,i])
    mu[4,i] <- (1-phi[3,i]) * (1-phi[2,i]) * (1-phi[1,i])
    for ( r in 1:(Nout-1) ) {
      phi[r,i] <- ilogit( zbeta0[r] + sum( zbeta[r,1:Nx] * zx[i,1:Nx] ) )
    }
  }
  for ( r in 1:(Nout-1) ) {
    zbeta0[r] ~ dnorm( 0 , 1/20^2 )
    for ( j in 1:Nx ) {
      zbeta[r,j] ~ dnorm( 0 , 1/20^2 )
    }
  }
  for ( r in 1:(Nout-1) ) {
    beta[r,1:Nx] <- zbeta[r,1:Nx] / xsd[1:Nx]
    beta0[r] <- zbeta0[r] - sum( zbeta[r,1:Nx] * xm[1:Nx] / xsd[1:Nx] )
  }
}
"
writeLines( model , "model.txt" )

dat<-read.csv( "http://www.matsuka.info/data_folder/SoftmaxRegData2.csv" )
y = dat$Y
x = cbind(dat[,1],dat[,2])
Ntotal = length(y)
Nout = length(unique(y))
dataList = list(y = y, x = x, Nx = 2, Ntotal = Ntotal, Nout = Nout)
parameters = c( "beta0" ,  "beta")
jagsModel = jags.model( "model.txt", data=dataList, n.chains=3, n.adapt=500 )
update( jagsModel , n.iter=1000)
codaSamples = coda.samples( jagsModel , variable.names=parameters, n.iter=10000, thin=5)
mcmcMat<-as.matrix(codaSamples)
par(mfrow=c(1,1))
plot(x[,1],x[,2],col=y)
n2plot=100
idx=sample(1:nrow(mcmcMat),n2plot)
temp.x = seq(-3,3,length.out = 100)
for (i.cat in 0:2){
  for (i_sample in 1:n2plot) {
    abline(a=-1*mcmcMat[idx[i_sample],7+i.cat]/mcmcMat[idx[i_sample],4+i.cat],
           b=-1*mcmcMat[idx[i_sample],1+i.cat]/mcmcMat[idx[i_sample],4+i.cat],col="orange")
  }
}
model2 = "
data {
  for ( j in 1:Nx ) {
    xm[j]  <- mean(x[,j])
    xsd[j] <-   sd(x[,j])
      for ( i in 1:Ntotal ) {
        zx[i,j] <- ( x[i,j] - xm[j] ) / xsd[j]
      }
    }
  }
model {
  for ( i in 1:Ntotal ) {
    y[i] ~ dcat( mu[1:Nout,i] )
    mu[1,i] <- phi[2,i] * phi[1,i]
    mu[2,i] <- (1-phi[2,i]) * phi[1,i]
    mu[3,i] <- phi[3,i] * (1-phi[1,i])
    mu[4,i] <- (1-phi[3,i]) * (1-phi[1,i])
    for ( r in 1:(Nout-1) ) {
      phi[r,i] <- ilogit( zbeta0[r] + sum( zbeta[r,1:Nx] * zx[i,1:Nx] ) )
    }
  }
  for ( r in 1:(Nout-1) ) {
    zbeta0[r] ~ dnorm( 0 , 1/20^2 )
    for ( j in 1:Nx ) {
      zbeta[r,j] ~ dnorm( 0 , 1/20^2 )
    }
  }
  for ( r in 1:(Nout-1) ) {
    beta[r,1:Nx] <- zbeta[r,1:Nx] / xsd[1:Nx]
    beta0[r] <- zbeta0[r] - sum( zbeta[r,1:Nx] * xm[1:Nx] / xsd[1:Nx] )
  }
}"
writeLines( modelString , con="TEMPmodel.txt" )

院 認識情報解析

library(rjags)
source("http://peach.l.chiba-u.ac.jp/course_folder/HDI_revised.txt")
dat = read.csv("http://peach.l.chiba-u.ac.jp/course_folder/HtWtData30.csv")

x = dat$height
y = dat$weight

modelS.txt = "
data{
  Ntotal <- length(y)
  xm <- mean(x)
  ym <- mean(y)
  xsd <- sd(x)
  ysd <- sd(y)
  for (i in 1:length(y)){
    zx[i] <- (x[i] - xm)/xsd
    zy[i] <- (y[i] - ym)/ysd
  }
}
model{
  for (i in 1:Ntotal){
    zy[i] ~ dt(zbeta0 + zbeta1 * zx[i], 1/zsigma^2, nu)
  }
  zbeta0 ~ dnorm(0, 1/10^2)
  zbeta1 ~ dnorm(0, 1/10^2)
  zsigma ~ dunif(1.03E-3, 1.03E+3)
  nu <- nuMinusOne + 1
  nuMinusOne ~ dexp(1/29.0)

  #Transfrom to original scale:
  beta1 <- zbeta1 * ysd/xsd
  beta0 <- zbeta0 * ysd + ym -zbeta1*xm*ysd/xsd
  sigma <- zsigma * ysd
}
"
writeLines(modelS.txt, "modelS.txt")
dataList = list(x = x, y = y)
jagsModel =jags.model("modelS.txt", data = dataList, n.chains = 3, n.adapt = 500)
update(jagsModel, 500)
codaSamples = coda.samples(jagsModel, variable.names = c("beta1", "beta0", "sigma","zbeta0","zbeta1"), n.iter = ((10000*1)/1), n.adapt = 500)
mcmcMat<-as.matrix(codaSamples)

plot(dat$height,dat$weight,xlab="height",ylab ='weight',pch=20,col="orange",cex=5)
par(mfrow=c(2,2))
HDI.plot(mcmcMat[,1])
plot(mcmcMat[,2],mcmcMat[,1], xlab='B1',ylab="B0",pch=20, col='orange')
plot(mcmcMat[,1],mcmcMat[,2], xlab='B0',ylab="B1",pch=20, col='orange')
HDI.plot(mcmcMat[,2])

par(mfrow=c(2,2))
HDI.plot(mcmcMat[,4])
plot(mcmcMat[,5],mcmcMat[,4], xlab='B1',ylab="B0",pch=20, col='orange')
plot(mcmcMat[,4],mcmcMat[,5], xlab='zB0',ylab="zB1",pch=20, col='orange')
HDI.plot(mcmcMat[,5])

n.mcmc = nrow(mcmcMat)
par(mfrow=c(1,1))
temp.x = c(0,80)
temp.y = mcmcMat[1,1]+mcmcMat[1,2]*temp.x
plot(temp.x,temp.y,type='l',lwd=2,col="orange",xlab="height",ylab='weight')
idx = sample(1:n.mcmc,100)
for (i.plot in 1:length(idx)){
  temp.y = mcmcMat[idx[i.plot],1]+mcmcMat[idx[i.plot],2]*temp.x
  lines(temp.x,temp.y,lwd=2,col="orange")
}
points(dat$height,dat$weight,pch=20,cex=4)


par(mfrow=c(1,1))
n2plot = 100
temp.x = c(0,80)
mean.set = seq(50,80,5)
idx=sample(1:nrow(mcmcMat),n2plot)
plot(x,y,xlim=c(50,80),ylim=c(0,400))
for (i_sample in 1:n2plot) {
  temp.y = mcmcMat[idx[i_sample],1]+mcmcMat[idx[i_sample],2]*temp.x
  lines(temp.x,temp.y,lwd=2,col="orange")
  for (i.means in 1:length(mean.set)){
    means = mcmcMat[idx[i_sample],1]+mcmcMat[idx[i_sample],2]*mean.set[i.means]
    y.seq = seq(0,400,length.out=1000)
    dens.y=dt((y.seq-means)/mcmcMat[idx[i_sample],3],29)
    dens.y=dens.y/max(dens.y)
    lines(dens.y,y.seq,type='l',col="orange")
    lines(-dens.y+mean.set[i.means],y.seq,type='l',col="orange")
  }
}
points(x,y,pch=20,cex=3)


model.txt = "
data{
  Ntotal <- length(y)
  xm <- mean(x)
  ym <- mean(y)
  xsd <- sd(x)
  ysd <- sd(y)
  for (i in 1:length(y)){
    zx[i] <- (x[i] - xm) / xsd
    zy[i] <- (y[i] - ym) / ysd
  }
}
model{
  for (i in 1:Ntotal){
    zy[i] ~ dt( zbeta0[s[i]] + zbeta1[s[i]] * zx[i], 1 / zsigma^2, nu)
  }
  for (j in 1:Nsubj){
    zbeta0[j] ~ dnorm( zbeta0mu, 1/(zbeta0sigma)^2)
    zbeta1[j] ~ dnorm( zbeta1mu, 1/(zbeta1sigma)^2)
  }

  zbeta0mu ~ dnorm(0, 1/(10)^2)
  zbeta1mu ~ dnorm(0, 1/(10)^2)
  zsigma ~ dunif(1.0E-3, 1.0E+3)
  zbeta0sigma ~ dunif(1.0E-3, 1.0E+3)
  zbeta1sigma ~ dunif(1.0E-3, 1.0E+3)
  nu <- nuMinusOne + 1
  nuMinusOne ~ dexp(1/29.0)

  for (j in 1:Nsubj){
    beta1[j] <- zbeta1[j] * ysd / xsd
    beta0[j] <- zbeta0[j] * ysd + ym - zbeta1[j] * xm * ysd / xsd
  }
  beta1mu <- zbeta1mu * ysd / xsd
  beta0mu <- zbeta0mu * ysd + ym -zbeta1mu * xm * ysd / xsd
  sigma <- zsigma * ysd
}
"
writeLines(model.txt, "model.txt")

dat = read.csv( file="http://peach.l.chiba-u.ac.jp/course_folder/HierLinRegressData.csv" )
x = dat$X
y = dat$Y
s = dat$Subj

dataList = list(x = x ,  y = y ,s = s , Nsubj = max(s))

jagsModel =jags.model("model.txt", data = dataList, n.chains = 3, n.adapt = 500)
update(jagsModel, 500)
codaSamples = coda.samples(jagsModel, variable.names = c("beta1mu", "beta0mu", "sigma",  "beta0","beta1"),
                           n.iter = ((10000*1)/1), n.adapt = 500)
mcmcMat<-as.matrix(codaSamples)

par(mfrow=c(1,1))

temp.x = c(40,90)
temp.y = mcmcMat[1,"beta0mu"]+mcmcMat[1,"beta1mu"]*temp.x
plot(temp.x,temp.y,type='l',lwd=2,col="orange",xlab="height",ylab='weight',xlim=c(40,90),ylim=c(40,260))
idx = sample(1:nrow(mcmcMat),100)
for (i.plot in 1:length(idx)){
  temp.y = mcmcMat[idx[i.plot],"beta0mu"]+mcmcMat[idx[i.plot],"beta1mu"]*temp.x
  lines(temp.x,temp.y,lwd=2,col="orange")
}
points(dat$X,dat$Y,pch=20,cex=4)


par(mfrow=c(5,5))
temp.x = c(40,90)
for (i.s in 1:25){
  temp.y = mcmcMat[1,i.s]+mcmcMat[1,i.s+26]*temp.x
  plot(temp.x,temp.y,type='l',lwd=2,col="orange",xlab="height",ylab='weight',xlim=c(40,90),ylim=c(40,260))
  idx = sample(1:nrow(mcmcMat),100)
  for (i.plot in 1:length(idx)){
    temp.y = mcmcMat[idx[i.plot],i.s]+mcmcMat[idx[i.plot],i.s+26]*temp.x
    lines(temp.x,temp.y,lwd=2,col="orange")
  }
  points(dat$X[which(dat$Subj==i.s)],dat$Y[which(dat$Subj==i.s)],pch=20,cex=2)
}

# ch 18
dat = read.csv( file="httP://www.matsuka.info/univ/course_folder/Guber1999data.csv" )
par(mfrow=c(1,1))
plot(dat[,c(2,5,8)],pch=20,cex=3)
y = dat$SATT
x = as.matrix(cbind(dat$Spend,dat$PrcntTake))


dataList = list(x = x ,y = y, Nx = dim(x)[2], Ntotal = dim(x)[1])

model.txt = "
data{
  ym <- mean(y)
  ysd <- sd(y)
  for (i in 1:Ntotal){
    zy[i] <- (y[i] - ym) / ysd
  }
  for (j in 1:Nx){ 
    xm[j] <- mean(x[,j]) 
    xsd[j] <- sd(x[,j])
    for (i in 1:Ntotal){
      zx[i,j] <- (x[i,j] - xm[j])/xsd[j]
    }
  }
}
model{
  for( i in 1:Ntotal){
    zy[i] ~ dt( zbeta0 + sum( zbeta[1:Nx] * zx[i, 1:Nx]), 1/zsigma^2, nu)
  }

  zbeta0 ~ dnorm(0, 1/2^2)
  for ( j in 1:Nx){
    zbeta[j] ~ dnorm(0, 1/2^2)
  }
  zsigma ~ dunif(1.0E-5, 1.0E+1)
  nu <- nuMinusOne + 1
  nuMinusOne ~ dexp(1/29.0)

  beta[1:Nx] <- ( zbeta[1:Nx] / xsd[1:Nx] ) * ysd
  beta0 <- zbeta0 * ysd + ym - sum( zbeta[1:Nx] * xm[1:Nx] / xsd[1:Nx]) * ysd
  sigma <- zsigma * ysd
}
"
writeLines(model.txt, "model.txt")

jagsModel =jags.model("model.txt", data = dataList, n.chains = 3, n.adapt = 500)
update(jagsModel, 500)
codaSamples = coda.samples(jagsModel, variable.names = c("beta[1]", "beta[2]", "beta0", "sigma"), n.iter = ((10000*1)/1), n.adapt = 500)
mcmcMat<-as.matrix(codaSamples)

plot(as.data.frame(mcmcMat))
par(mfrow=c(1,2))
HDI.plot(mcmcMat[,2])
HDI.plot(mcmcMat[,3])

x = as.matrix(cbind(dat$Spend,dat$PrcntTake,dat$Spend*dat$PrcntTake))

dataList = list(x = x ,y = y, Nx = dim(x)[2], Ntotal = dim(x)[1])
jagsModel =jags.model("model.txt", data = dataList, n.chains = 3, n.adapt = 500)
update(jagsModel, 500)
codaSamples = coda.samples(jagsModel, variable.names = c("beta[1]", "beta[2]", "beta[3]", "beta0", "sigma"), n.iter = ((10000*1)/1), n.adapt = 500)
mcmcMat<-as.matrix(codaSamples)
par(mfrow=c(1,3))
HDI.plot(mcmcMat[,2])
HDI.plot(mcmcMat[,3])
HDI.plot(mcmcMat[,4])

empHDI(mcmcMat[,2])

par(mfrow=c(1,1))
perT = seq(0,80,5)
plot(x=c(perT[1],perT[1]),empHDI(mcmcMat[,2]+mcmcMat[,4]*perT[1]),
     xlim=c(-5,85),ylim=c(-15,35),type='l',lwd=3,col="orange",
     xlab = c("value on Percent Take"),ylab="Slope on Spend")
for (i.plot in 2:length(perT)){
  lines(c(perT[i.plot],perT[i.plot]),empHDI(mcmcMat[,2]+mcmcMat[,4]*perT[i.plot]),
        lwd=3,col="orange")
}
abline(h=0,lwd=2,lty=2)

par(mfrow=c(1,1))
perT = seq(0,10,0.5)
plot(x=c(perT[1],perT[1]),empHDI(mcmcMat[,3]+mcmcMat[,4]*perT[1]),
     xlim=c(-0.5,11),ylim=c(-5,1),type='l',lwd=3,col="orange",
     xlab = c("value on Spent"),ylab="Slope on %Take")
for (i.plot in 2:length(perT)){
  lines(c(perT[i.plot],perT[i.plot]),empHDI(mcmcMat[,3]+mcmcMat[,4]*perT[i.plot]),
        lwd=3,col="orange")
}
abline(h=0,lwd=2,lty=2)

院:認識情報解析 T-test & ANOVA

source("http://www.matsuka.info/univ/course_folder/HDI_revised.txt")
dat <- read.csv("http://peach.l.chiba-u.ac.jp/course_folder/TwoGroupIQ.csv")
dat2sd <- dat[dat$Group=="Smart Drug",]

dataList <- list(y = dat2sd$Score, Ntotal = nrow(dat2sd),
                 meanY = mean(dat2sd$Score), sdY = sd(dat2sd$Score))

model.txt = "
model {
  for (i in 1:Ntotal){
    y[i] ~ dnorm(mu, 1/sigma^2)
  }
  mu ~ dnorm(meanY, 1/(100*sdY)^2)
  sigma ~ dunif(sdY/1000, sdY*1000)
}"
writeLines(model.txt, "model.txt")
parameters = c( "mu" ,  "sigma")
library("rjags")
jagsModel = jags.model( "model.txt", data=dataList, n.chains=3, n.adapt=500 )
update( jagsModel , n.iter=1000)
codaSamples = coda.samples( jagsModel , variable.names=parameters, n.iter=10000, thin=5)
traceplot(codaSamples)
gelman.plot(codaSamples)
mcmcMat<-as.matrix(codaSamples)

HDI.plot(mcmcMat[,1])
HDI.plot(mcmcMat[,2])
hist(dat2sd$Score, probability = T)
x.temp = seq(40,220,0.1)
n.plot = 100
rand.order = sample(1:nrow(mcmcMat),n.plot)
for (i.plot in 1:n.plot){
  y.temp = dnorm(x.temp, mean = mcmcMat[rand.order[i.plot],1],
                 sd = mcmcMat[rand.order[i.plot],2])
  lines(x.temp, y.temp, type='l',col="orange")
}

#y = c(-2:2,15);Ntotal = length(y);meanY = mean(y); sdY = sd(y)
#dataList = list(y= y, Ntotal=Ntotal, meanY = meanY, sdY = sdY)
dataList <- list(y = dat2sd$Score, Ntotal = nrow(dat2sd),
                 meanY = mean(dat2sd$Score), sdY = sd(dat2sd$Score))

model.txt ="
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 <- nuMinusOne+1
  nuMinusOne ~ dexp(1/29)
}"
writeLines(model.txt, "model.txt")
parameters = c( "mu" ,  "sigma", "nu")
jagsModel = jags.model( "model.txt", data=dataList, n.chains=3, n.adapt=1000 )
update( jagsModel , n.iter=1000)
codaSamples = coda.samples( jagsModel , variable.names=parameters, n.iter=10000, thin=5)
#traceplot(codaSamples)
#gelman.plot(codaSamples)
mcmcMat<-as.matrix(codaSamples)
HDI.plot(mcmcMat[,1])
HDI.plot(mcmcMat[,2])

dataList <- list(y = dat$Score, Ntotal = nrow(dat),
                 meanY = mean(dat$Score), sdY = sd(dat$Score),
                 Ngroup = length(unique(dat$Group)),
                 group = dat$Group)

model.txt ="
model {
  for ( i in 1:Ntotal ) {
    y[i] ~ dt( mu[group[i]] , 1/sigma[group[i]]^2 , nu )
  }
  for (j in 1:Ngroup){
    mu[j] ~ dnorm( meanY , 1/(100*sdY)^2 )
    sigma[j] ~ dunif( sdY/1000 , sdY*1000 )
  }
  nu <- nuMinusOne+1
  nuMinusOne ~ dexp(1/29)
}"
writeLines(model.txt, "model.txt")
parameters = c( "mu" ,  "sigma", "nu")
jagsModel = jags.model( "model.txt", data=dataList, n.chains=3, n.adapt=1000 )
update( jagsModel , n.iter=1000)
codaSamples = coda.samples( jagsModel , variable.names=parameters, n.iter=10000, thin=5)
#traceplot(codaSamples)
#gelman.plot(codaSamples)
mcmcMat<-as.matrix(codaSamples)
HDI.plot(mcmcMat[,1]-mcmcMat[,2])
HDI.plot(mcmcMat[,2])


dat<-read.csv("http://www.matsuka.info/univ/course_folder/FruitflyDataReduced.csv")
y=dat$Longevity
yMean = mean(y)
ySD = sd(y)
Ntotal = length(y)

MeanSD2gamma <- function( mean, sd ) {
  shape = mean^2 / sd^2
  rate = mean / sd^2
  return(data.frame(shape,rate))
}

ModeSD2gamma <- function( mode, sd ) {
  rate = ( mode + sqrt( mode^2 + 4 * sd^2 ) )/( 2 * sd^2 )
  shape = 1 + mode * rate
  return(data.frame(shape,rate))
}

temp.gamma = ModeSD2gamma( mode=sd(y)/2 , sd=2*sd(y) )
gShape = temp.gamma$shape
gRate = temp.gamma$rate

x=as.numeric(dat$CompanionNumber)
Ngroup = length(unique(x))
xc=dat$Thorax
xcMean=mean(xc)

dataList = list(
  y = y ,
  x = x ,
  Ntotal = Ntotal ,
  Ngroup = Ngroup ,
  yMean = yMean ,
  ySD = ySD ,
  gShape = gShape,
  gRate  = gRate
)

model.txt="
model {
  for ( i_data in 1:Ntotal ) {
    y[ i_data ] ~ dnorm( a0 + a[ x[ i_data ] ], 1/ySigma^2)
  }

  ySigma ~ dunif( ySD/100, ySD*10 )
  a0 ~ dnorm( yMean, 1/(ySD*5)^2)

  for ( i_group in 1:Ngroup ) {
    a[ i_group ] ~ dnorm(0.0, 1/aSigma^2)
  }
  aSigma ~ dgamma( gShape, gRate )

  for ( i_group in 1:Ngroup ) { m[ i_group ] <- a0 + a[ i_group ] }
  b0 <- mean( m[ 1:Ngroup ] )
  for (i_data in 1:Ngroup ) { b[ i_data ] <- m[ i_data ] - b0 }
  }
"
writeLines(model.txt, "model.txt")

parameters = c( "b0" ,  "b" ,  "ySigma",  "aSigma" )

jagsModel = jags.model( "model.txt", data=dataList, n.chains=3, n.adapt=500 )
update( jagsModel , n.iter=1000)
codaSamples = coda.samples( jagsModel , variable.names=parameters, n.iter=10000, thin=5)
mcmcMat<-as.matrix(codaSamples)

# fig 19.3
plot(x,dat$Longevity,xlim=c(0,5.1),ylim=c(0,120))
for (i_c in 1:100) {
  sd=mcmcMat[i_c,"ySigma"]
  lim=qnorm(c(0.025,0.975))
  means=mcmcMat[i_c,2:6]+mcmcMat[i_c,"b0"]
  for (i_group in 1:5){
    lower=means[i_group]+lim[1]*sd
    upper=means[i_group]+lim[2]*sd
    yseq=seq(lower,upper,length.out = 100)
    dens.y=dnorm((yseq-means[i_group])/sd)
    dens.y=dens.y/max(dens.y)*0.75
    lines(-dens.y+i_group,yseq,type='l',col="orange")
  }
}




dataList = list(
  y = y ,
  x = x ,
  xc = xc ,
  Ntotal = Ntotal ,
  Ngroup = Ngroup ,
  yMean = yMean ,
  ySD = ySD ,
  xcMean = xcMean ,
  acSD = sd(xc) ,
  gShape = gShape,
  gRate  = gRate
)
model.txt="
model {
  for ( i_data in 1:Ntotal ) {
    y[ i_data ] ~ dnorm( mu[i_data], 1/ySigma^2)
    mu [ i_data] <- a0 + a[ x[ i_data ] ] + ac*( xc[ i_data ] - xcMean )
  }

  ySigma ~ dunif( ySD/100, ySD*10 )
  a0 ~ dnorm( yMean, 1/(ySD*5)^2)

  for ( i_group in 1:Ngroup ) {
    a[ i_group ] ~ dnorm(0.0, 1/aSigma^2)
  }
  aSigma ~ dgamma( gShape, gRate )
  ac ~ dnorm(0, 1/(2*ySD/acSD)^2 )
  b0 <- a0 + mean( a[ 1:Ngroup ] ) - ac*xcMean
  for (i_data in 1:Ngroup ) { b[ i_data ] <- a[ i_data ] - mean( a[1:Ngroup]) }
}"
writeLines(model.txt, "model.txt")
parameters = c( "b0" ,  "b" ,  "ySigma",  "aSigma", "ac" )

jagsModel = jags.model( "model.txt", data=dataList, n.chains=3, n.adapt=500 )
update( jagsModel , n.iter=5000)
codaSamples = coda.samples( jagsModel , variable.names=parameters, n.iter=10000, thin=5)
mcmcMat<-as.matrix(codaSamples)

# fig 19.5 only "none0" will be plotted
plot(Longevity~Thorax,xlim=c(0.6,1),ylim=c(0,120),data=dat[dat$CompanionNumber=="None0",],type='n')
Thorax.value=c(0.65,0.8,0.95)
for (i_c in 1:100) {

  # regression lines
  xs=c(0.5,1.1)
  ys=mcmcMat[i_c,"b0"]+mcmcMat[i_c,3]+mcmcMat[i_c,"ac"]*xs
  lines(xs,ys,col="orange")

  # density
  sd=mcmcMat[i_c,"ySigma"]
  lim=qnorm(c(0.025,0.975))
  means=mcmcMat[i_c,"b0"]+mcmcMat[i_c,3]+mcmcMat[i_c,"ac"]*Thorax.value
  for (i_thorax in 1:3){
    lower=means[i_thorax]+lim[1]*sd
    upper=means[i_thorax]+lim[2]*sd
    yseq=seq(lower,upper,length.out = 100)
    dens.y=dnorm((yseq-means[i_thorax])/sd)
    dens.y=dens.y/max(dens.y)*0.05
    lines(-dens.y+Thorax.value[i_thorax],yseq,type='l',col="orange")
  }
}
abline(v=Thorax.value,col='green')
points(Longevity~Thorax,xlim=c(0.6,1),ylim=c(0,120),data=dat[dat$CompanionNumber=="None0",])



dat<-read.csv("http://www.matsuka.info/univ/course_folder/NonhomogVarData.csv")
y=dat$Y
yMean = mean(y)
ySD = sd(y)
Ntotal = length(y)
x=as.numeric(dat$Group)
Ngroup = length(unique(x))

model.txt = "
model {
  for ( i_data in 1:Ntotal ) {
y[ i_data ] ~ dt( a0 + a[ x[ i_data ] ], 1/ySigma[x[i_data]]^2, nu)
}
nu <- nuMinusOne + 1
nuMinusOne ~ dexp(1/29)
for (i_group in 1:Ngroup) {
ySigma[i_group]~dgamma(ySigSh,ySigR)
}
ySigSh <- 1+ySigMode * ySigR
ySigR  <- ((ySigMode + sqrt(ySigMode^2+4*ySigSD^2))/(2*ySigSD^2))

ySigMode ~ dgamma(gShape,gRate)
ySigSD ~ dgamma(gShape,gRate)
a0 ~ dnorm( yMean, 1/(ySD*10)^2)

for ( i_group in 1:Ngroup ) {
a[ i_group ] ~ dnorm(0.0, 1/aSigma^2)
}
aSigma ~ dgamma( gShape, gRate )

for ( i_group in 1:Ngroup ) {
m[ i_group ] <- a0 + a[ i_group ]
}

b0 <- mean( m[ 1:Ngroup ] )

for (i_data in 1:Ngroup ) {
b[ i_data ] <- m[ i_data ] - b0
}
}"
writeLines(model.txt, "model.txt")

dataList = list(
  y = y ,
  x = x ,
  Ntotal = Ntotal ,
  Ngroup = Ngroup ,
  yMean = yMean ,
  ySD = ySD ,
  gShape = gShape,
  gRate  = gRate
)

parameters = c( "b0" ,  "b" ,  "ySigma",  "aSigma" ,"nu","ySigMode","ySigSD")

jagsModel = jags.model( "model.txt", data=dataList, n.chains=3, n.adapt=500 )
update( jagsModel , n.iter=1000)
codaSamples = coda.samples( jagsModel , variable.names=parameters, n.iter=10000, thin=5)
mcmcMat<-as.matrix(codaSamples)


# plotting
plot(x,y,xlim=c(0,4.1),ylim=c(70,130))
n2plot=100
idx=sample(1:nrow(mcmcMat),n2plot)
for (i_sample in 1:n2plot) {
  i_c = idx[i_sample]
  lim=qnorm(c(0.025,0.975))
  means=mcmcMat[i_c,2:6]+mcmcMat[i_c,"b0"]
  for (i_group in 1:4){
    sd=mcmcMat[i_c,9+i_group]
    lower=means[i_group]+lim[1]*sd
    upper=means[i_group]+lim[2]*sd
    yseq=seq(lower,upper,length.out = 100)
    dens.y=dnorm((yseq-means[i_group])/sd)
    dens.y=dens.y/max(dens.y)*0.75
    lines(-dens.y+i_group,yseq,type='l',col="orange")
  }
}

院:認知情報解析

y = sample(c(rep(1,15), rep(0,35)))
Ntotal=length(y)
datalist = list(y=y,Ntotal=Ntotal)
source("http://www.matsuka.info/univ/course_folder/HDI_revised.txt")

library(rjags)
txt = "
model {
  for ( i_data in 1:Ntotal ) {
    y[ i_data ] ~ dbern( theta )
  }
  theta ~ dbeta( 1, 1 )
}"
writeLines(txt, "~/model.txt")
jagsModel = jags.model(file="~/model.txt",
                       data=datalist,n.chains=3,n.adapt=500)
update(jagsModel,n.iter=1000)
codaSamples=coda.samples(jagsModel,variable.names=c("theta"),n.iter=5000)
mcmcMat<-as.matrix(codaSamples)

HDI.plot(mcmcMat)
traceplot(codaSamples)
autocorr.plot(codaSamples,type='l')
gelman.plot(codaSamples)

y1 = sample(c(rep(1,6), rep(0,2)))
y2 = sample(c(rep(1,2), rep(0,5))) 
y = c(y1,y2)
s = c(rep(1,8),rep(2,7))
Ntotal=length(y)
datalist = list(y = y, Ntotal = Ntotal, s = s, Nsubj = 2)

txt = "
model {
for ( i_data in 1:Ntotal ) {
y[ i_data ] ~ dbern( theta[s[i_data]] )
}
for ( i_s in 1:Nsubj) {
theta[i_s] ~ dbeta( 2, 2 )
}
}"
writeLines(txt, "~/model.txt")
jagsModel = jags.model(file="~/model.txt",
                       data=datalist,n.chains=3,n.adapt=500)
update(jagsModel,n.iter=1000)
codaSamples=coda.samples(jagsModel,variable.names=c("theta"),n.iter=5000)
mcmcMat<-as.matrix(codaSamples)
# checking MCMC
HDI.plot(mcmcMat[,2])
traceplot(codaSamples)
autocorr.plot(codaSamples)
gelman.plot(codaSamples)

x.temp = seq(0,150,0.1)
g1 = dgamma(x.temp, shape =0.01, rate =0.01)
plot(x.temp,g1,type='l')

g2 = dgamma(x.temp, shape =1.56, rate = 0.0312)
plot(x.temp,g2,type='l')

g3 = dgamma(x.temp, shape =6.25, rate = 0.125)
plot(x.temp,g2,type='l')

txt = "
model {
for ( i_data in 1:Ntotal ) {
y[ i_data ] ~ dbern( theta[s[i_data]] )
}
for ( i_s in 1:Nsubj) {
theta[i_s] ~ dbeta( omega*(kappa-2)+1,(1-omega)*(kappa-2)+1)
}
omega ~ dbeta(1,1)
kappa <- kappaMinusTwo + 2
kappaMinusTwo ~ dgamma(0.01, 0.01)
}"
writeLines(txt, "~/model.txt")


dat<-read.csv("http://www.matsuka.info/data_folder/TherapeuticTouchData.csv")
y=dat$y
s=as.numeric(dat$s)
Ntotal=length(dat$y)
Nsubj=length(unique(s))
datalist = list(y=y,s=s,Ntotal=Ntotal,Nsubj=Nsubj)

jagsModel = jags.model(file="~/model.txt",data=datalist,n.chains=3,n.adapt=500)
update(jagsModel,n.iter=1000)
codaSamples=coda.samples(jagsModel,variable.names=c("theta","omega","kappa"),n.iter=5000)
mcmcMat<-as.matrix(codaSamples)

dat<-read.csv("http://www.matsuka.info/data_folder/BattingAverage.csv")
z = dat$Hits
N = dat$AtBats
s = dat$PlayerNumber
c = dat$PriPosNumber
Ncat = 9
Nsubj = 948

datalist = list(z = z, N=N, s=s, c=c, Ncat=Ncat, Nsubj =Nsubj )
txt = "
model {
 for (i_s in 1:Nsubj) {
   z[i_s] ~ dbin( theta[i_s], N[i_s])
   theta[i_s] ~ dbeta(omega[c[i_s]]*(kappa[c[i_s]]-2)+1,
      (1-omega[c[i_s]])*(kappa[c[i_s]]-2)+1)
 }
 for (i_c in 1:Ncat) {
   omega[i_c] ~ dbeta(omegaO*(kappaO-2)+1,(1-omegaO)*(kappaO-2)+1)
   kappa[i_c] <- kappaMinusTwo[i_c]+2
   kappaMinusTwo[i_c] ~ dgamma(0.01,0.01)
 }
 omegaO ~ dbeta(1,1)
 kappaO <- kappaMinusTwoO+2
 kappaMinusTwoO ~ dgamma(0.01,0.01)
}"
writeLines(txt, "~/model.txt")
jagsModel = jags.model(file="~/model.txt",data=datalist,n.chains=3,n.adapt=500)
update(jagsModel,n.iter=1000)
codaSamples=coda.samples(jagsModel,variable.names=c("theta","omega","kappa"),n.iter=5000)
mcmcMat<-as.matrix(codaSamples)
HDI.plot(mcmcMat[," omega[1]"]) #pitcher
HDI.plot(mcmcMat[," omega[2]"]) #catcher
HDI.plot(mcmcMat[," omega[1]"]) #1st base

院:認知情報解析学

source("http://www.matsuka.info/univ/course_folder/cuUtil02.R")
dat<-read.csv('http://www.matsuka.info/data_folder/tdkCFA.csv')
dat.fa1 <- factanal(dat,1)
dat.fa2 <- factanal(dat,2)
dat.fa3 <- factanal(dat,3)
dat.fa4 <- factanal(dat,4)

install.packages('sem')
library(sem)

model01=cfa(reference.indicator=FALSE)
F1:extrovert,cheerful, leadership, antisocial, talkative, motivated, hesitance, popularity

mod1<-sem(model01, data = dat)

model02=cfa(reference.indicator=FALSE)
F1: extrovert, leadership, motivated, hesitance 
F2: cheerful, antisocial, talkative, popularity

mod2<-sem(model02, data = dat)

認知情報解析学演習 RNN

txt = "You say goodbye and I say hello.  you say goodbay and I say hello"

txt2corpus <- function(txt){
  txt = tolower(txt)
  txt = gsub('[.]', ' . sos',txt)
  words = unlist(strsplit(c('sos',txt), " "))
  uniq.words = unique(words)
  n.uniq = length(uniq.words)
  n.words = length(words)
  corpus = rep(0,n.words)
  corpus = match(words,uniq.words)
  return(corpus)
}

corp = txt2corpus(txt)

corp2contxt1SRNN = function(corpus){
  len.corp = length(corpus)
  # creating target matrix
  idxT = cbind(1:(len.corp-1), corpus[2:len.corp])
  targ1S = matrix(0,nrow=len.corp-1,ncol=length(unique(corpus)))
  targ1S[idxT]=1
  # creating context matrices
  idxC = cbind(1:(len.corp-1),corpus[1:(len.corp-1)])
  contxt = matrix(0,nrow=len.corp-1,ncol=length(unique(corpus)))
  contxt[idxC]=1
  return(list(y=targ1S,x=contxt))
}

init.RNN <- function(n.uniq,size.hidden){
  W.h = matrix(rnorm(size.hidden*size.hidden),nrow=size.hidden)*0.01
  W.x = matrix(rnorm(n.uniq*size.hidden),nrow=n.uniq)*0.01
  b = matrix(rnorm(size.hidden),nrow=1)*0.01
  return(list(W.h = W.h, W.x= W.x, b = b))
}
MatMult.forwd <- function(x, W){
  return(x%*%W)
}

MatMult.bckwd <- function(x, W, dout){
  dx = dout%*%t(W)
  dW = t(x)%*%dout
  return(list(dx = dx, dW = dW))
}

affine.forwd <- function(x, W, b){
  return(x%*%W + matrix(1, nrow = nrow(x), ncol = 1)%*%b)
}

affine.bckwd <- function(x, W, b, dout){
  dx = dout%*%t(W)
  dW = t(x)%*%dout
  db = colSums(dout)
  return(list(dx = dx, dW = dW, db = db))
}


softmax.forwd <- function(x, target){
  max.x = apply(x,1,max)
  C = ncol(x)
  x = x - max.x%*%matrix(1,nrow=1,ncol=C)
  y = exp(x)/rowSums(exp(x))
  delta = 1e-7;
  R = nrow(as.matrix(y))
  return(-sum(target*log(y + delta))/R)
}

softmax.pred <- function(x, target){
  max.x = apply(x,1,max)
  C = ncol(x)
  x = x - max.x%*%matrix(1,nrow=1,ncol=C)
  y = exp(x)/rowSums(exp(x))
  return(y)
}

softmax.bckwd <- function(x, target,  dout = 1){
  max.x = apply(x, 1, max)
  R = nrow(x)
  C = ncol(x)
  x = x - max.x%*%matrix(1,nrow=1,ncol=C)
  y = exp(x)/rowSums(exp(x))
  return((y-target)/R)
}


RNN.forward <- function(h.prev, x, network){
  b.size = nrow(x)
  h = h.prev%*%network$W.h + x%*%network$W.x
  hb = h+matrix(1,nrow=b.size,ncol=1)%*%network$b
  h.next = tanh(hb)
  return(h.next = h.next)
}

RNN.backward <- function(dh.next, network, x, h.next, h.prev){
  dt = dh.next * (1- h.next^2)
  db = colSums(dt)
  dW.h = t(h.prev)%*%dt
  dh.prev = dt%*%t(network$W.h)
  dW.x = t(x)%*%dt
  dx = dt%*%t(network$W.x)
  return(list(db = db, dW.h = dW.h, dh.prev = dh.prev, dW.x = dW.x, dx=dx))
}


txt = "You say goodbye and I say hello."

batch.size = 3
time = 3

size.hidden = 7
corp = txt2corpus(txt)
dat<-corp2contxt1SRNN(corp)
network <- init.RNN(ncol(dat$y), size.hidden)
corp.len = nrow(dat$y)
h.prev =array(0, c(batch.size, size.hidden, time))
h.next = array(0, c(batch.size, size.hidden, (time+1)))
W.out = matrix(rnorm(size.hidden * ncol(dat$y)), nrow = size.hidden)
b.out = matrix(rnorm(ncol(dat$y)), nrow = 1)
n.rep = 100000;lambda = 0.25;

loss = rep(0,n.rep)
for (i.rep in 1:n.rep){
  for (i.t in 1:time){
    idx = seq(i.t, corp.len, time)
    h.next[, , i.t] = RNN.forward(h.prev[, , i.t], dat$x[idx,], network)
    if (i.t < time){
      h.prev[, , (i.t+1)] = h.next[, , i.t]
    }
  }
  dWHs = matrix(0,nrow=nrow(network$W.h),ncol=ncol(network$W.h))
  dWXs = matrix(0,nrow=nrow(network$W.x),ncol=ncol(network$W.x))
  dBs = matrix(0,nrow=nrow(network$b),ncol=ncol(network$b))
  dWOs = matrix(0,nrow=nrow(W.out),ncol=ncol(W.out))
  dBOs = matrix(0,nrow=nrow(b.out),ncol=ncol(b.out))
  d.prev = matrix(0,nrow=batch.size,ncol=size.hidden)
  L = 0
  for (i.t in time:1){
    idx = idx = seq(i.t, corp.len, time)
    O = affine.forwd(h.next[,,i.t], W.out, b.out)
    L = L + softmax.forwd(O, dat$y[idx,])
    ds = softmax.bckwd(O, dat$y[idx,], 1)
    dW.o = affine.bckwd(h.next[,,i.t], W.out, b.out, ds)
    dWOs = dWOs + dW.o$dW
    dBOs = dBOs + dW.o$db
    RNN.d = RNN.backward(dW.o$dx+d.prev, network, dat$x[idx,],h.next[,,(i.t+1)],h.prev[,,i.t])
    #RNN.d = RNN.backward(dW.o$dx, network, dat$x[idx,],h.next[,,i.t],h.prev[,,i.t])
    dWHs = dWHs + RNN.d$dW.h
    dWXs = dWXs + RNN.d$dW.x
    dBs = dBs + RNN.d$db
    d.prev = RNN.d$dh.prev
  }
  loss[i.rep] = L
  W.out = W.out - lambda*dWOs
  b.out = b.out - lambda*dBOs
  network$W.h - lambda*dWHs
  network$W.x - lambda*dWXs
  network$b - lambda*dBs
}
plot(loss, type='l')
for (i.t in 1:time){
  idx = idx = seq(i.t, corp.len, time)
  O = affine.forwd(h.next[,,i.t], W.out, b.out)
  print(softmax.pred(O, dat$y[idx,]))
}

院:認知情報解析学演習 ch23 & 24

dat<-read.csv("http://peach.l.chiba-u.ac.jp/course_folder/OrdinalProbitData-1grp-1.csv" )

model = "
model {
  for ( i in 1:Ntotal ) {
    y[i] ~ dcat( pr[i,1:nYlevels] )
    pr[i,1] <- pnorm( thresh[1] , mu , 1/sigma^2 )
    for ( k in 2:(nYlevels-1) ) {
      pr[i,k] <- max( 0 ,  pnorm( thresh[ k ] , mu , 1/sigma^2 )- pnorm( thresh[k-1] , mu , 1/sigma^2 ) )
    }
    pr[i,nYlevels] <- 1 - pnorm( thresh[nYlevels-1] , mu , 1/sigma^2 )
  }
  mu ~ dnorm( (1+nYlevels)/2 , 1/(nYlevels)^2 )
  sigma ~ dunif( nYlevels/1000 , nYlevels*10 )
  for ( k in 2:(nYlevels-2) ) {  # 1 and nYlevels-1 are fixed, not stochastic
    thresh[k] ~ dnorm( k+0.5 , 1/2^2 )
  }
}
"
writeLines( model , "model.txt" )

y = dat$Y
Ntotal = length(y)
nYlevels = max(y)
thresh = rep(NA,nYlevels-1)
thresh[1] = 1 + 0.5
thresh[nYlevels-1] = nYlevels-1 + 0.5
dataList = list(y = y , nYlevels = nYlevels,thresh = thresh, Ntotal = Ntotal)

parameters = c( "pr" ,  "thresh","mu","sigma")
jagsModel = jags.model( "model.txt", data=dataList, n.chains=3, n.adapt=500 )
update( jagsModel , n.iter=1000)
codaSamples = coda.samples( jagsModel , variable.names=parameters, n.iter=10000, thin=5)
mcmcMat<-as.matrix(codaSamples)

meanT = rowMeans(cbind(mcmcMat[,"thresh[1]"],mcmcMat[,"thresh[2]"],mcmcMat[,"thresh[3]"],mcmcMat[,"thresh[4]"],mcmcMat[,"thresh[5]"],
                   mcmcMat[,"thresh[6]"]))
plot(mcmcMat[,"thresh[1]"],meanT,xlim=c(0.5,7))
points(mcmcMat[,"thresh[2]"],meanT,xlim=c(0.5,7))
points(mcmcMat[,"thresh[3]"],meanT,xlim=c(0.5,7))
points(mcmcMat[,"thresh[4]"],meanT,xlim=c(0.5,7))
points(mcmcMat[,"thresh[5]"],meanT,xlim=c(0.5,7))
points(mcmcMat[,"thresh[6]"],meanT,xlim=c(0.5,7))
HDI.plot(mcmcMat[,"mu"],xlab="mu")
HDI.plot(mcmcMat[,"sigma"],xlab="sigma")


model2 = "
  model {
for ( i in 1:Ntotal ) {
y[i] ~ dcat( pr[i,1:nYlevels] )
pr[i,1] <- pnorm( thresh[1] , mu[x[i]] , 1/sigma[x[i]]^2 )
for ( k in 2:(nYlevels-1) ) {
pr[i,k] <- max( 0 ,  pnorm( thresh[ k ] , mu[x[i]] , 1/sigma[x[i]]^2 )
- pnorm( thresh[k-1] , mu[x[i]] , 1/sigma[x[i]]^2 ) )
}
pr[i,nYlevels] <- 1 - pnorm( thresh[nYlevels-1] , mu[x[i]] , 1/sigma[x[i]]^2 )
}
for ( j in 1:2 ) { # 2 groups
mu[j] ~ dnorm( (1+nYlevels)/2 , 1/(nYlevels)^2 )
sigma[j] ~ dunif( nYlevels/1000 , nYlevels*10 )
}
for ( k in 2:(nYlevels-2) ) {  # 1 and nYlevels-1 are fixed, not stochastic
thresh[k] ~ dnorm( k+0.5 , 1/2^2 )
}
}
"
writeLines( model2 , "model2.txt" )

dat<-read.csv("http://peach.l.chiba-u.ac.jp/course_folder/OrdinalProbitData1.csv")
y = dat$Y
x = as.numeric(dat$X)
Ntotal = length(y)
nYlevels = max(y)
thresh = rep(NA,nYlevels-1)
thresh[1] = 1 + 0.5
thresh[nYlevels-1] = nYlevels-1 + 0.5
dataList = list(y = y , x= x ,nYlevels = nYlevels,thresh = thresh, Ntotal = Ntotal)

parameters = c( "pr" ,  "thresh","mu","sigma")
jagsModel = jags.model( "model2.txt", data=dataList, n.chains=3, n.adapt=500 )
update( jagsModel , n.iter=1000)
codaSamples = coda.samples( jagsModel , variable.names=parameters, n.iter=10000, thin=5)
mcmcMat<-as.matrix(codaSamples)
HDI.plot(mcmcMat[,"mu[1]"])
HDI.plot(mcmcMat[,"mu[2]"])
HDI.plot(mcmcMat[,"mu[1]"]-mcmcMat[,"mu[2]"])


model3 = "
data {
xm  <- mean(x)
xsd <-   sd(x)
for ( i in 1:Ntotal ) {
zx[i] <- ( x[i] - xm ) / xsd
}
}
model {
for ( i in 1:Ntotal ) {
y[i] ~ dcat( pr[i,1:nYlevels] )
pr[i,1] <- pnorm( thresh[1] , mu[i] , 1/sigma^2 )
for ( k in 2:(nYlevels-1) ) {
pr[i,k] <- max( 0 ,  pnorm( thresh[ k ] , mu[i] , 1/sigma^2 )
- pnorm( thresh[k-1] , mu[i] , 1/sigma^2 ) )
}
pr[i,nYlevels] <- 1 - pnorm( thresh[nYlevels-1] , mu[i] , 1/sigma^2 )
mu[i] <- zbeta0 + sum( zbeta * zx[i] )
}
# Priors vague on standardized scale:
zbeta0 ~ dnorm( (1+nYlevels)/2 , 1/(nYlevels)^2 )
zbeta ~ dnorm( 0 , 1/(nYlevels)^2 )

zsigma ~ dunif( nYlevels/1000 , nYlevels*10 )
# Transform to original scale:
beta <- ( zbeta / xsd )
beta0 <- zbeta0  - sum( zbeta * xm / xsd )
sigma <- zsigma
for ( k in 2:(nYlevels-2) ) {  # 1 and nYlevels-1 are fixed
thresh[k] ~ dnorm( k+0.5 , 1/2^2 )
}
}
"
writeLines( model3, con="model3.txt" )

dat<-read.csv("http://peach.l.chiba-u.ac.jp/course_folder/OrdinalProbitData-LinReg-2.csv")
y = dat$Y
x = as.numeric(dat$X)
Ntotal = length(y)
Nx =1
nYlevels = max(y)
thresh = rep(NA,nYlevels-1)
thresh[1] = 1 + 0.5
thresh[nYlevels-1] = nYlevels-1 + 0.5
dataList = list(y = y , Nx= 1,x= x ,nYlevels = nYlevels,thresh = thresh, Ntotal = Ntotal)

parameters = c("thresh","mu","sigma","beta","beta0")
jagsModel = jags.model( "model3.txt", data=dataList, n.chains=3, n.adapt=500 )
update( jagsModel , n.iter=1000)
codaSamples = coda.samples( jagsModel , variable.names=parameters, n.iter=10000, thin=5)
mcmcMat<-as.matrix(codaSamples)
plot(dat$X,dat$Y,pch=20,cex=3)
for (i.plot in 1:100){
  abline(a=mcmcMat[i.plot,"beta0"],b=mcmcMat[i.plot,"beta"],col="orange",lwd=2)
}
HDI.plot(mcmcMat[,"beta"])
#plot(1,1,type='n',xlim=c(0,6),ylim=c(0,6))
#abline(a=0,b=1,lwd=2)
#x.temp = seq(0,6,length.out=100)
#y = dnorm(x.temp,mean=3,sd=1)
#lines(x=-y*3+1,y=x.temp)
#lines(x=c(-1,3),y=c(3,3))
#lines(x=c(-1,0.5),y=c(0.5,0.5),lty=2,col="red")
#lines(x=c(-1,1.5),y=c(1.5,1.5),lty=2,col="red")
#lines(x=c(-1,2.5),y=c(2.5,2.5),lty=2,col="red")
#lines(x=c(-1,3.5),y=c(3.5,3.5),lty=2,col="red")
#lines(x=c(-1,4.5),y=c(4.5,4.5),lty=2,col="red")
#lines(x=c(-1,5.5),y=c(5.5,5.5),lty=2,col="red")

model4="
model {
  for ( i in 1:Ncell ) {
    y[i] ~ dpois( lambda[i] )
    lambda[i] <- exp( a0 + a1[x1[i]] + a2[x2[i]] + a1a2[x1[i],x2[i]] )
  }
  a0 ~ dnorm( yLogMean , 1/(yLogSD*2)^2 )
  for ( j1 in 1:Nx1Lvl ) { a1[j1] ~ dnorm( 0.0 , 1/a1SD^2 ) }
  a1SD ~ dgamma(agammaShRa[1],agammaShRa[2])
  for ( j2 in 1:Nx2Lvl ) { a2[j2] ~ dnorm( 0.0 , 1/a2SD^2 ) }
  a2SD ~ dgamma(agammaShRa[1],agammaShRa[2])
  for ( j1 in 1:Nx1Lvl ) { for ( j2 in 1:Nx2Lvl ) {
    a1a2[j1,j2] ~ dnorm( 0.0 , 1/a1a2SD^2 )
  } }
  a1a2SD ~ dgamma(agammaShRa[1],agammaShRa[2])
  # Convert a0,a1[],a2[],a1a2[,] to sum-to-zero b0,b1[],b2[],b1b2[,] :
  for ( j1 in 1:Nx1Lvl ) { for ( j2 in 1:Nx2Lvl ) {
    m[j1,j2] <- a0 + a1[j1] + a2[j2] + a1a2[j1,j2] # cell means
  } }
  b0 <- mean( m[1:Nx1Lvl,1:Nx2Lvl] )
  for ( j1 in 1:Nx1Lvl ) { b1[j1] <- mean( m[j1,1:Nx2Lvl] ) - b0 }
  for ( j2 in 1:Nx2Lvl ) { b2[j2] <- mean( m[1:Nx1Lvl,j2] ) - b0 }
  for ( j1 in 1:Nx1Lvl ) { for ( j2 in 1:Nx2Lvl ) {
    b1b2[j1,j2] <- m[j1,j2] - ( b0 + b1[j1] + b2[j2] )
  } }
  # Compute predicted proportions:
  for ( j1 in 1:Nx1Lvl ) { for ( j2 in 1:Nx2Lvl ) {
    expm[j1,j2] <- exp(m[j1,j2])
    ppx1x2p[j1,j2] <- expm[j1,j2]/sum(expm[1:Nx1Lvl,1:Nx2Lvl])
  } }
  for ( j1 in 1:Nx1Lvl ) { ppx1p[j1] <- sum(ppx1x2p[j1,1:Nx2Lvl]) }
  for ( j2 in 1:Nx2Lvl ) { ppx2p[j2] <- sum(ppx1x2p[1:Nx1Lvl,j2]) }
}"
writeLines( model4, con="model4.txt" )

dat = read.csv( file="http://peach.l.chiba-u.ac.jp/course_folder/HairEyeColor.csv" )

y=dat$Count
x1=dat$Eye
x2=dat$Hair

x1 = as.numeric(as.factor(x1))
x1levels = levels(as.factor(dat$Eye))
x2 = as.numeric(as.factor(x2))
x2levels =  levels(as.factor(dat$Eye))
Nx1Lvl = length(unique(x1))
Nx2Lvl = length(unique(x2))
Ncell = length(y) # number of rows of data, not sum(y)

yLogMean = log(sum(y)/(Nx1Lvl*Nx2Lvl))
yLogSD = log(sd(c(rep(0,Ncell-1),sum(y))))
MeanSD2gamma <- function( mean, sd ) {
  shape = mean^2 / sd^2
  rate = mean / sd^2
  return(data.frame(shape,rate))
}

ModeSD2gamma <- function( mode, sd ) {
  rate = ( mode + sqrt( mode^2 + 4 * sd^2 ) )/( 2 * sd^2 )
  shape = 1 + mode * rate
  return(data.frame(shape,rate))
}
temp=ModeSD2gamma(mode=yLogSD , sd=2*yLogSD)
agammaShRa = unlist(temp )

data_list = list(
  y = y ,x1 = x1 , x2 = x2 ,
  Ncell = Ncell , Nx1Lvl = Nx1Lvl , Nx2Lvl = Nx2Lvl ,
  yLogMean = yLogMean , yLogSD = yLogSD ,agammaShRa = agammaShRa
)

jagsModel =jags.model("model4.txt", data = data_list, n.chains = 3, n.adapt = 500)
update(jagsModel, 500)
codaSamples = coda.samples(jagsModel, variable.names = c("b0", "b1", "b2", "b1b2", "ppx1p", "ppx2p", "ppx1x2p"),
                           n.iter = ((10000*1)/1), n.adapt = 500)
mcmcMat<-as.matrix(codaSamples)
EyeBlueHairBlack <- mcmcMat[,"ppx1x2p[1,1]"]
HDI.plot(EyeBlueHairBlack)


EyeBlue_vs_EyeBrown_at_HairBlack <- mcmcMat[,"b1b2[1,1]"] - mcmcMat[,"b1b2[2,1]"] - mcmcMat[,"b1[2]"] + mcmcMat[,"b1[1]"]
HDI.plot(EyeBlue_vs_EyeBrown_at_HairBlack)

EyeBlue_vs_EyeBrown_at_Blond <- mcmcMat[,"b1b2[1,2]"] - mcmcMat[,"b1b2[2,2]"] - mcmcMat[,"b1[2]"] + mcmcMat[,"b1[2]"]
HDI.plot(EyeBlue_vs_EyeBrown_at_Blond)

diff <- EyeBlue_vs_EyeBrown_at_HairBlack - EyeBlue_vs_EyeBrown_at_Blond
HDI.plot(diff)

認知情報解析演習

# data preparation
dat<-read.csv("~/Downloads/DBDA2Eprograms/z15N50.csv")
y=dat$y
Ntotal=length(dat$y)
datalist = list(y=y,Ntotal=Ntotal)

# model
model.txt = "
model {
  for ( i_data in 1:Ntotal ) {
    y[ i_data ] ~ dbern( theta )
  }
  theta ~ dbeta( 1, 1 )
}
"
writeLines(model.txt, "model.txt")

# jags
library(rjags)
jagsModel = jags.model(file="model.txt",data=datalist,n.chains=3,n.adapt=500)
update(jagsModel,n.iter=1000)
codaSamples=coda.samples(jagsModel,variable.names=c("theta"),n.iter=5000)
mcmcMat<-as.matrix(codaSamples)

par(mfrow=c(2,2))
cols=c("orange", "skyblue","pink")

# chain
par(mfrow=c(2,2))
cols=c("orange", "skyblue","pink")
mcmc1<-as.mcmc(codaSamples[[1]])
mcmc2<-as.mcmc(codaSamples[[2]])
mcmc3<-as.mcmc(codaSamples[[3]])
plot(as.matrix(mcmc1),type='l',col=cols[1])
lines(as.matrix(mcmc2),col=cols[2])
lines(as.matrix(mcmc3),,col=cols[3])

# autocorrelation
ac1=autocorr(mcmc1,lags=0:50)
ac2=autocorr(mcmc2,lags=0:50)
ac3=autocorr(mcmc3,lags=0:50)
plot(ac1, type='o', pch = 20, col = cols[1], ylab = "Autocorrelation", xlab = "Lag")
lines(ac2, type='o', pch = 20, col = cols[2])
lines(ac3, type='o', pch = 20, col = cols[3])

# shrink factor
resALL=mcmc.list(mcmc1,mcmc2,mcmc3)
gd1=gelman.plot(resALL, auto.layout = F)

# density
den1=density(mcmc1)
den2=density(mcmc2)
den3=density(mcmc3)
plot(den1, type='l', col = cols[1], main = " ", xlab = "param. value",lwd=2.5)
lines(den2, col = cols[2], lwd=2.5)
lines(den3, col = cols[3], lwd=2.5)

par(mfrow=c(1,1))

認知情報解析


### model 
model {
  for ( i in 1:Ntotal ) {
    y[i] ~ dnorm( mu , 1/sigma^2  )
  }
  mu ~ dnorm( meanY , 1/(100*sdY)^2 )
  sigma ~ dunif( sdY/1000 , sdY*1000 )
}
###

dat<-read.csv("~/Downloads/DBDA2Eprograms/TwoGroupIQ.csv")
y = dat$Score[dat$Group=="Smart Drug"]
Ntotal = length(y)
dataList = list(y = y, Ntotal = Ntotal, meanY = mean(y), sdY = sd(y))


dataList = list(y = y, Ntotal = Ntotal, meanY = mean(y), sdY = sd(y))
jagsModel = jags.model("~/research/oka/DBDA/ch16/model_16_1_2.txt", data=dataList, n.chains=3, n.adapt=1000 )
update( jagsModel , n.iter=1000)
codaSamples = coda.samples( jagsModel , variable.names=c("mu","sigma"), n.iter=5000, thin=5)
mcmcMat<-as.matrix(codaSamples)

hist(y,nclass=15,probability = T)
x.temp = seq(40,200,0.1)
n.plot = 100
randperm = sample(1:nrow(mcmcMat),n.plot)
for (i.plot in 1:n.plot){
  norm.temp=dnorm(x.temp,mean=mcmcMat[randperm[i.plot],1],sd=mcmcMat[randperm[i.plot],2])
  lines(x.temp,norm.temp,col='orange')
}

### model 16.2
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 <- nuMinusOne+1
  nuMinusOne ~ dexp(1/30.0)
}
###
hist(y,nclass=20,probability = T)
n.plot = 100
randperm = sample(1:nrow(mcmcMat),n.plot)
for (i.plot in 1:n.plot){
  x.temp1 = seq(40,200,0.1)
  x.temp2 = (x.temp1 - mcmcMat[randperm[i.plot],1])/mcmcMat[randperm[i.plot],3]
  t.temp=dt(x.temp2,mcmcMat[randperm[i.plot],2])/mcmcMat[randperm[i.plot],3]
  lines(x.temp1,t.temp,col='orange')
}

### model 16.3
model {
  for ( i in 1:Ntotal ) {
    y[i] ~ dt( mu[group[i]] , 1/sigma[group[i]]^2 , nu )
  }
  for (j in 1:Ngroup){
    mu[j] ~ dnorm( meanY , 1/(100*sdY)^2 )
    sigma[j] ~ dunif( sdY/1000 , sdY*1000 )
  }      
  nu <- nuMinusOne+1
  nuMinusOne ~ dexp(1/29)
}
###
dat<-read.csv("~/Downloads/DBDA2Eprograms/TwoGroupIQ.csv")
y = dat$Score
group = as.numeric(dat$Group)
Ntotal = length(y)
Ngroup = length(unique(group))

dataList = list(y = y, Ntotal = Ntotal, meanY = mean(y), sdY = sd(y),Ngroup=Ngroup,group=group)
jagsModel = jags.model("~/research/oka/DBDA/ch16/model_16_3.txt", data=dataList, n.chains=3, n.adapt=5000 )
update( jagsModel , n.iter=1000)
codaSamples = coda.samples( jagsModel , variable.names=c("mu","sigma","nu"), n.iter=5000, thin=5)
mcmcMat<-as.matrix(codaSamples)