x=rnorm(n=1, mean=100, sd=15)
y=runif(n=3, min=1, max=10)
N = 10000
random.data = rnorm(N, mean=0, sd=1)
hist(random.data, nclass = 50, col = "navy", xlab = "Data",
probability = T, main = "Histogram of Random Data")
dens = density(random.data)
lines(dens, col = "orange", lwd = 4)
sample(1:10,3)
sample(c(“gu”,“choki”,“pa”),1)
sample(1:10)
sample(0:1, 10, replace=T)
# FOR loop
for (i_loop in 1:5) {
print(i_loop)
}
for (i_loop in 1:5) {
print(c(i_loop, 2^i_loop))
}
counter <- 1
while(counter<=10){
print(counter)
counter<-counter+1
}
counter <- 1
while(counter^2 <= 10){
print(c(counter, counter^2))
counter<-counter+1
}
affil<-"cogsci"
if (affil=="cogsci") {
print("you are wonderful")
}
affil<-"phil"
if (affil=="cogsci") {
print("you are wonderful")
} else {
print("still, you are wonderful")
}
counter=6
repeat{
print(counter)
counter = counter + 1
if(counter>5){break}
}
counter=6
repeat{
if(counter>5){break}
print(counter)
counter+counter+1
}
six.counter=0; N = 1000
for (i_loop in 1:N) {
die<-sample(1:6,1)
if (die==6) {six.counter=six.counter+1}
}
six.counter/N
N=1000; six.counter=rep(0,N);
for (i_loop in 1:N) {
die<-sample(1:6,1)
if (die==6) {six.counter[i_loop]=1}
}
plot(1:N,cumsum(six.counter)/(1:1000),type='l',ylim=c(0,1),
lwd=2)
abline(h=1/6,lwd=2,col='red')
N = 1000
die.all <- sample(1:6,N,replace=T)
six.index <- die.all==6
par(mfrow = c(2,1))
par(oma=c(2,2,0,0),mar=c(4,4,1,1),mfrow=c(2,1))
plot(1:N, die.all, pch=20, col = 'red', ylim = c(0,7),
ylab = "Result", xlab = "trial")
plot(1:N,cumsum(six.index)/(1:1000), type='l', ylim=c(0,1), lwd=2,
ylab = "P(die = 6)", xlab = "trial")
abline(h=1/6,lwd=2,col='red')
# CLT w/ loop
N=10;nRep=10000;
means<-rep(0,nRep)
for (i_rep in 1:nRep) {
dat<-runif(N)
means[i_rep]=mean(dat)
}
hist(means,nclass=50,probability=T)
dens<-density(means)
lines(dens,col='skyblue',lwd=3)
xs=seq(-0,1,0.01)
theo.dens<-dnorm(xs,mean=0.5,sd=sqrt((1/12)/N))
lines(xs,theo.dens,col='orange',lwd=3,lty=2)
# CLT w/o loop
N=10
nRep=10000
dat<-matrix(runif(N*nRep),nrow=N)
means<-colMeans(dat)
hist(means,nclass=50,probability=T)
dens<-density(means); lines(dens,col='skyblue',lwd=3)
xs=seq(-0,1,0.01)
theo.dens<-dnorm(xs,mean=0.5,sd=sqrt((1/12)/N))
lines(xs,theo.dens,col='orange',lwd=3,lty=2)
r=-99;v1=1633;v2=355
while (r!=0){
r=v1%%v2
print(paste('v1 =',v1,', v2 = ',v2,',remainder = ',r))
v1=v2
v2=r
}
tol = 1e-7;grad = 1e10; lambda = 0.1;x = 10; x.hist = x
grad = 2*x+2
while (abs(grad)>tol){
x = x - lambda*grad
x.hist=c(x.hist,x)
grad = 2*x+2
}
x.temp=seq(-10,10,0.1)
plot(x.temp, x.temp^2+2*x.temp,type='l',lwd=2)
lines(x.hist,x.hist^2+2*x.hist,type='o',pch=1,col='red',cex=1)
データ解析基礎論a 分散分析３
source("http://www.matsuka.info/univ/course_folder/cutil.R")
dat<-read.csv("http://www.matsuka.info/data_folder/dktb316.txt")
dat.aov<-aov(words~method+subj+Error(subj:method),data=dat)
dat.aov2<-aov(words~method+Error(subj+subj:method),data=dat)
dat.aov.BTW <-aov(words~method,data=dat)
summary(dat.aov)
RB.qtukey(dat.aov,dat, 0.05)
dat<-read.csv("http://www.matsuka.info/data_folder/dktb3211.txt")
interaction.plot(dat$duration, # x軸
dat$method, # まとめる変数
dat$result, # y軸
pch=c(20,20),
col=c("skyblue","orange"),
ylab="score",
xlab="Duration",
lwd=3,
type='b',
cex=2,
trace.label="Method")
dat.aov <- aov(result~method*duration+Error(s+s:duration),dat)
TSME<-SPF.tsme(dat.aov,dat,"result")
dat<-read.csv("http://www.matsuka.info/data_folder/dktb3218.txt")
dat.aov<-aov(result~method+duration+method:duration + Error(s+method:s+duration:s+method:duration:s), data=dat)
TSME<-RBF.tsme(dat.aov, dat, "result")
認知情報解析演習 居住区問題
n.circle = 20; n.sharp = 20; size = 10
loc = sample(1:size^2, n.circle+n.sharp)
type = c(rep(1,n.circle),rep(2,n.sharp))
# circle = 1; sharp = 2
people = cbind(type,loc)
state.mat = matrix(0,size,size)
state.mat[people[,2]]=people[,1]
p.move = #1/(2*n.circle)
p.other = 0.1
c.count = 3
idx = cbind(rep(1:size,size),sort(rep(1:size,size)))
#
term = F
while (term == F){
term = T
rand.order = sample(1:nrow(people))
for (i.p in 1:nrow(people)){
if (people[rand.order[i.p],1]==1){
# circle
if (runif(1) < p.move){
empty = 1:(size^2)
empty = empty[-sort(people[,2])]
people[rand.order[i.p],2] = sample(empty,1)
state.mat = matrix(0,size,size)
state.mat[people[,2]]=people[,1]
term = F
}
} else {
# sharp
x.min = max(idx[people[rand.order[i.p],2],1]-1,1)
x.max = min(idx[people[rand.order[i.p],2],1]+1,size)
y.min = max(idx[people[rand.order[i.p],2],2]-1,1)
y.max = min(idx[people[rand.order[i.p],2],2]+1,size)
circle.in = sum(state.mat[x.min:x.max,y.min:y.max]==1)
if (circle.in >= c.count){
empty = 1:(size^2)
empty = empty[-sort(people[,2])]
people[rand.order[i.p],2] = sample(empty,1)
state.mat = matrix(0,size,size)
state.mat[people[,2]]=people[,1]
term = F
#print('moved')
}
}
}
for (i.p in 1:nrow(people)){
if (people[rand.order[i.p],1]==2){
x.min = max(idx[people[i.p,2],1]-1,1)
x.max = min(idx[people[i.p,2],1]+1,size)
y.min = max(idx[people[i.p,2],2]-1,1)
y.max = min(idx[people[i.p,2],2]+1,size)
circle.in = sum(state.mat[x.min:x.max,y.min:y.max]==1)
print(circle.in)
if (circle.in >= c.count){
term = F
break
}
}
}
}
plot(0,0, type= 'n', xlim = c(0,11),ylim=c(0,11))
lab = c("0","#")
text(idx[people[,2],1],idx[people[,2],2],lab[people[,1]],cex=3)
ab = seq(0.5,10.5,1)
for (i in 1:11){
abline(h=ab[i],col='red')
abline(v=ab[i],col='red')
}
認知情報解析演習 石とりゲーム（bugあり）
col1 = matrix(c(rep(0,4),c(1,0,0,0),c(1,1,0,0),c(1,1,1,0),rep(1,4)),nrow=4,byrow=T)
col2 = matrix(c(rep(10,4),c(11,10,10,10),c(11,11,10,10),c(11,11,11,10),rep(11,4)),nrow=4,byrow=T)
col3 = matrix(c(rep(100,4),c(101,100,100,100),c(101,101,100,100),c(101,101,101,100),rep(101,4)),nrow=4,byrow=T)
act.list = list()
state.temp = list()
counter = 0
Q1 = list()
Q2 = list()
for (i.c1 in 1:5){
if (sum(col1[,i.c1])==0){
act1 = c()
} else {
act1 = seq(1,sum(col1[,i.c1]),1)
}
for (i.c2 in 1:5){
if (sum(col2[,i.c2])==40){
act2 = c()
} else {
act2 = seq(11,sum(col2[,i.c2]==11)*11,11)
}
for (i.c3 in 1:5){
if (sum(col3[,i.c3])==400){
act3 = c()
} else {
act3 = seq(101,sum(col3[,i.c3]==101)*101,101)
}
counter = counter + 1
state.temp[[counter]] = cbind(col1[,i.c1],col2[,i.c2],col3[,i.c3])
act.list[[counter]] = c(act1,act2,act3)
Q1[[counter]] = rep(0, length(c(act1,act2,act3)))
Q2[[counter]] = rep(0, length(c(act1,act2,act3)))
}
}
}
rm.stone <- function(act, st.shape){
if (act == -99){s}
if (act > 100){
n.remove = act%%100
n.stone = length(which(st.shape[,3]==101))
start = (5-n.stone)
st.shape[(start:(start+n.remove-1)),3] = 100
} else {
if (act > 10){
n.remove = act%%10
n.stone = length(which(st.shape[,2]==11))
start = (5-n.stone)
st.shape[(start:(start+n.remove-1)),2] = 10
} else {
n.remove = act
n.stone = length(which(st.shape[,1]==1))
start = (5-n.stone)
st.shape[(start:(start+n.remove-1)),1] = 0
}
}
return(st.shape)
}
id.state <- function(st.shape, state.temp){
for (i.st in 1:125){
if (all(st.shape == state.temp[[i.st]])){
state.idx = i.st
break
}
}
return(state.idx)
}
ck.act <- function(Q, act.vec, eta){
if (is.null(act.vec)){
return(list(act = -99, act.idx = -99))
break
}
if (length(act.vec)==1){
act = act.vec
} else {
p = exp(Q[[state.idx]])/sum(exp(Q[[state.idx]]))
act = sample(act.vec, 1, prob = p)
}
act.idx = which(act.vec==act)
return(list(act = act, act.idx = act.idx))
}
gamma=1;alpha = 0.1;n.rep=10000
for (i.rep in 1:n.rep){
# first action
state.idx = 125; counter = 1
st.shape = state.temp[[state.idx]]
res.act = ck.act(Q1,act.list[[state.idx]],eta)
act = res.act$act;act.idx = res.act$act.idx
state.old = state.idx
act.old = act.idx
# 2nd to last
while (state.idx != 1) {
counter = counter + 1
st.shape <- rm.stone(act, st.shape)
state.idx <- id.state(st.shape, state.temp)
if (counter%%2==1) {
res.act = ck.act(Q1,act.list[[state.idx]],eta)
} else {
res.act = ck.act(Q2,act.list[[state.idx]],eta)
}
act = res.act$act; act.idx = res.act$act.idx
if (state.idx == 1){
if (counter%%2==1){rew1 = -1; rew2 = 1;} else {rew1 = 1; rew2 = -1;}
Q1[[state.old]][act.old]=Q1[[state.old]][act.old]
+alpha*(rew1-Q1[[state.old]][act.old])
Q2[[state.old]][act.old]=Q2[[state.old]][act.old]
+alpha*rew2-Q2[[state.old]][act.old])
} else {
rew1 = 0;
rew2 =0;
if (counter%%2==1){
Q1[[state.old]][act.old]=Q1[[state.old]][act.old]
+alpha*(rew1+gamma* Q1[[state.idx]][act.idx]-Q1[[state.old]][act.old])
} else {
Q2[[state.old]][act.old]=Q2[[state.old]][act.old]
+alpha*(rew2+gamma* Q2[[state.idx]][act.idx]-Q2[[state.old]][act.old])
}
}
state.old = state.idx
act.old = act.idx
}
}
2019年度 データ解析基礎論a 分散分析２
source("http://www.matsuka.info/univ/course_folder/cutil.R")
adj.alpha(5,0.05)
# anova
dat<-read.csv("http://www.matsuka.info/data_folder/dktb312.csv")
dat$method=factor(dat$method, levels(dat$method)[c(1,3,4,2)])
dat.aov<-aov(result~method, data=dat)
summary(dat.aov)
# multiple comparison
# do not use these command
# use simpler one - "cu.bonF1F"
dat.means<-tapply(dat$result,dat$method,mean)
new.alpha = 1-(1-0.05)^(1/5)
cont=c(-3,1,1,1)
bunshi=sum(cont*dat.means)
bunbo=sqrt(5.29*(sum((cont^2)/8)))
t.value=bunshi/bunbo
2*(1-pt(t.value,28))
# cu.bonF1F
new.alpha<-adj.alpha(5,0.05)
cu.bonF1F(dat.aov,dat,c(-3,1,1,1),new.alpha)
cu.bonF1F(dat.aov,dat,c(-1,1,0,0),new.alpha)
cu.bonF1F(dat.aov,dat,c(-1,0,1,0),new.alpha)
cu.bonF1F(dat.aov,dat,c(-1,0,0,1),new.alpha)
cu.bonF1F(dat.aov,dat,c(0,-2,1,1),new.alpha)
# cu.scheffe1F
new.f<-adj.f(4,3,28,0.05)
cu.scheffe1F(dat.aov,dat,c(-3,1,1,1))
# 2 way ANOVA
dat <- read.csv("http://peach.l.chiba-u.ac.jp/course_folder/datW03R.txt")
interaction.plot(dat$gender,
dat$affil,
dat$shoesize,
pch=c(20,20),
col=c("skyblue","orange"),
xlab="gender", ylab="shoesize",
lwd=3,type='b',cex=2,
trace.label="Affiliation")
dat.aov=aov(shoesize~gender*affil, data=dat)
dat.aov.sum=summary(dat.aov)
# testing simple main effect
# do not use these command
# use simpler one - CRF.tsme
means<-tapply(dat$shoesize, list(dat$gender,dat$affil), mean)
SS_gen_CS<- 5*(means[2,1]^2 + means[1,1]^2 -0.5*sum(means[,1])^2) # SS_gender CS
SS_gen_PS<- 5*(means[2,2]^2 + means[1,2]^2 -0.5*sum(means[,2])^2) # SS_gender PS
dat.aov.sum=summary(dat.aov) # ANOVA table
MSe=dat.aov.sum[[1]][4,3] # MSE from ANOVA table or MSe=0.62
dfE=dat.aov.sum[[1]][4,1] # DF for error or dfE=16
dfG=1 # DF for gender
F_gen_CS=(SS_gen_CS/dfG)/MSe # F-value for gender effect given CS
F_gen_PS=(SS_gen_PS/dfG)/MSe # F-value for gender effect given PS
P_gen_CS=1-pf(F_gen_CS,1,dfE) # p-value for gender effect given CS
P_gen_PS=1-pf(F_gen_PS,1,dfE)
SS_affil_F<- 5*(means[1,1]^2+means[1,2]^2-0.5*sum(means[1,])^2) #SS_affil | F
SS_affil_M<- 5*(means[2,1]^2+means[2,2]^2-0.5*sum(means[2,])^2) #SS_affil | M
dfA=1 # DF for affil
F_affil_F=SS_affil_F/dfA/MSe # F-value for affiliation effect | F
F_affil_M=SS_affil_M/dfA/MSe # F-value for affiliation effect | M
P_affil_F=1-pf(F_affil_F,1,dfE) # p-value for affiliation effect | F
P_affil_M=1-pf(F_affil_M,1,dfE) # p-value for affiliation effect | M
# testint simple main effect w/ CRF.tsme
tsme = CRF.tsme(dat.aov, dat)
# another 2-way ANOVA
dat <-read.csv("http://www.matsuka.info/data_folder/dktb321.txt")
interaction.plot(dat$duration,dat$method,dat$result, pch=c(20,20),
col=c("skyblue","orange"), ylab="score", xlab="Duration",
lwd=3,type='b',cex=2,trace.label="Method")
mod1=aov(result~method+duration,data=dat)
mod1.sum=print(summary(mod1))
mod2=aov(result~method*duration,data=dat)
mod2.sum=print(summary(mod2))
CRF.tsme(mod2, dat)
# plotting
dat<-read.csv("http://www.matsuka.info/univ/course_folder/datW03R.txt",header=T)
means<-tapply(dat$shoesize,list(dat$gender, dat$affil),mean)
Ns<-tapply(dat$shoesize,list(dat$gender, dat$affil),length)
sds<-tapply(dat$shoesize,list(dat$gender, dat$affil),sd)
sems<-sds/sqrt(Ns)
plot(c(0,1),means[,1],type='o',col='skyblue', xlim=c(-0.1,1.1), lwd=2, cex=2, pch=20,
ylim=c(min(means)*0.975, max(means)*1.025), xlab="gender", ylab="shoesize", xaxt="n")
axis(1,c(0,1),c("female","male"))
lines(c(0,1),means[,2],type='o',col='orange',lwd=2,cex=2,pch=20)
lines(c(0,0),c(means[1,1]-sems[1,1],means[1,1]+sems[1,1]),col="skyblue",lwd=2)
lines(c(0,0),c(means[1,2]-sems[1,2],means[1,2]+sems[1,2]),col="orange",lwd=2)
lines(c(1,1),c(means[2,1]-sems[2,1],means[2,1]+sems[2,1]),col="skyblue",lwd=2)
lines(c(1,1),c(means[2,2]-sems[2,2],means[2,2]+sems[2,2]),col="orange",lwd=2)
legend("topleft",c("CogSci","PsySci"),col=c("skyblue","orange"),lwd=2)
院：認識情報解析
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" )
2019 データ解析基礎論A 分散分析１
f=c(24.1,23.9,24.4,24.4,23.5)
m=c(25.6,26.1,25.8,25.9,26)
# ANOVA w/o aov
G.mean=mean(c(f,m))
ss.total=sum((c(f,m)-G.mean)^2)
ss.tr=sum((mean(f)-G.mean)^2)*5+sum((mean(m)-G.mean)^2)*5
ss.error=sum((f-mean(f))^2)+sum((m-mean(m))^2)
ss.tr+ss.error
df.tr = 2 - 1
df.error = (4-1)*2
ms.tr = ss.tr /df.tr
ms.error = ss.error /df.error
1-pf(f.value,1,8)
# ANOVA w/ aov
dat<-data.frame(ssize=c(f,m),gender=c(rep("f",5),rep("m",5)))
dat.aov<-aov(ssize ~ gender, data=dat)
summary(dat.aov)
dat<-read.table("http://www.matsuka.info/data_folder/aov01.txt")
dat.aov<-aov(shoesize~club, dat)
# TukeyHSD w/o TukeyHSD command
qT<-qtukey(0.95, 3, 67)
HSD<-qT*sqrt((2.243*(1/23+1/24))/2)
means<-tapply(dat$shoesize,dat$club,mean)
abs(outer(means,means,"-"))
abs(outer(means,means,"-"))>HSD
# TukeyHSD
TukeyHSD(dat.aov)
dat<-read.csv("http://www.matsuka.info/data_folder/dktb312.csv")
dat.aov<-aov(result~method, data=dat)
summary(dat.aov)
source("http://www.matsuka.info/univ/course_folder/cutil.R")
adj.alpha(5,0.05)
# multiple t-test
dat.means<-tapply(dat$result,dat$method,mean)
new.alpha = 1-(1-0.05)^(1/5)
cont=c(-3,1,1,1)
bunshi=sum(cont*dat.means)
bunbo=sqrt(5.29*(sum((cont^2)/8)))
t.value=bunshi/bunbo
2*(1-pt(t.value,28)) > new.alpha
# bonferroni w/ cUTIL
cu.bonF1F(dat.aov,dat,c(-3,1,1,1),new.alpha)
cu.bonF1F(dat.aov,dat,c(-1,1,0,0),new.alpha)
cu.bonF1F(dat.aov,dat,c(-1,0,1,0),new.alpha)
cu.bonF1F(dat.aov,dat,c(-1,0,0,1),new.alpha)
cu.bonF1F(dat.aov,dat,c(0,-2,1,1),new.alpha)
# scheffe w/ cUTIL
new.f<-adj.f(4,3,28,0.05)
cu.scheffe1F(dat.aov,dat,c(-3,1,1,1))
cu.scheffe1F(dat.aov,dat,c(-1,0,1,0))
cu.scheffe1F(dat.aov,dat,c(-1,0,0,1))
基礎実習 MA0２
select.act <- function(state){
poss.act = which(state=="N")
if (length(poss.act)==1){
act = poss.act
} else {
act = sample(poss.act, 1)
}
return(act)
}
ck.term <- function(state) {
term = F
result = "undecided"
term.idx = matrix(c(1,2,3,4,5,6,7,8,9,1,4,7,
2,5,8,3,6,9,1,5,9,3,5,7),ncol=3,byrow=T)
if (sum(state == "N")==0) {
term=T
result = "tie"
}
for (i.ck in 1:8) {
O.won = all(state[term.idx[i.ck,]]=="O")
X.won = all(state[term.idx[i.ck,]]=="X")
if (O.won ==1 ) {
term= T
result = "O won"
break
}
if (X.won ==1){
term=T
result = "X won"
break
}
}
return(list(term = term, result=result))
}
ck.term2 <- function(state) {
term = F
result = "undecided"
term.idx = matrix(c(1,2,3,4,5,6,7,8,9,1,4,7,
2,5,8,3,6,9,1,5,9,3,5,7),ncol=3,byrow=T)
if (sum(state == "N")==0) {
term=T
result = "tie"
}
O.won = apply(term.idx,1,function(x) all(state[x]=="O"))
X.won = apply(term.idx,1,function(x) all(state[x]=="X"))
if (any(O.won) == T ) {
term= T
result = "O won"
} else if (any(X.won) ==1){
term=T
result = "X won"
}
return(list(term = term, result=result))
}
tictactoe <-function(){
coord = matrix(c(1,3,1,2,1,1,2,3,2,2,2,1,3,3,
3,2,3,1), nrow = 9, byrow=T)
plot(0,0, type="n", xlim= c(0.5,3.5), ylim=c(0.5,3.5))
abline(h = 1.5); abline(h = 2.5)
abline(v = 1.5); abline(v = 2.5)
state = rep("N",9); marker = c("O","X")
repeat {
for (i.player in 1:2) {
act = select.act(state)
state[act] = marker[i.player]
text(coord[act,1],coord[act,2], marker[i.player], cex=4)
res = ck.term2(state)
if (res$term == T) { break }
Sys.sleep(0.5)
}
if (res$term == T) {
print(paste("result:", res$result))
break
}
}
}
# 実行例
> tictactoe()
[1] "result: X won"
# 対戦版 - 無作為に動くので非常に弱いです。
# 上のselect.actを変更することでマシになると思います。
play.tictactoe.R <-function(){
coord = matrix(c(1,3,1,2,1,1,2,3,2,2,2,1,3,3,
3,2,3,1), nrow = 9, byrow=T)
plot(coord, pch=paste(1:9), col='gray', cex=3, xlim= c(0.5,3.5), ylim=c(0.5,3.5))
abline(h = 1.5); abline(h = 2.5)
abline(v = 1.5); abline(v = 2.5)
state = rep("N",9);
repeat {
for (i.player in 1:2){
if (i.player == 1){
act = as.numeric(readline(prompt="Enter box ID: "))
} else {act = select.act(state)}
state[act] = marker[i.player]
text(coord[act,1],coord[act,2], marker[i.player], cex=4)
res = ck.term2(state)
if (res$term == T) { break }
Sys.sleep(0.5)
}
if (res$term == T) {
print(paste("result:", res$result))
break
}
}
}
# 実行例：
> play.tictactoe.R()
Enter box ID: 5
Enter box ID: 1
Enter box ID: 9
[1] "result: O won"
データ解析基礎論 回帰分析３
dat <- read.csv("http://www.matsuka.info/data_folder/hwsk8-17-6.csv")
dat.lm<-lm(ani~otouto, dat)
plot(hatvalues(dat.lm))
dat.lm$residuals
x = dat$otouto
mean.x = mean(x)
h = 1/length(x)+(x-mean.x)^2/sum((x - mean.x)^2)
hatvalues(dat.lm)
e.std = e/(slm$sigma*sqrt(1-h))
rstandard(dat.lm)
(e^2*h)/(2*slm$sigma^2*(1-h)^2)
cooks.distance(dat.lm)
dat<-read.csv("http://www.matsuka.info/data_folder/tdkReg02.csv")
plot(dat)
dat.lm<-lm(sales~., data=dat)
install.packages("DAAG")
library(DAAG)
vif(dat.lm)
lm.price<-lm(price~material+design+dump, data=dat)
p.rsq = summary(lm.price)$r.squared
vif.p = 1/(1-p.rsq)
summary(dat.lm)
dat<-read.table("http://www.matsuka.info/data_folder/tdkPATH01.txt",header=TRUE)
plot(dat,pch=20,cex=2)
dat.lm1<-lm(absence~interest,dat)
dat.lm2<-lm(study~interest+knowledge,dat)
dat.lm3<-lm(grade~knowledge+study+absence,dat)
院 認識情報解析
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)