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)
Category Archives: 院・情報解析
院：認識情報解析
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)