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)
Monthly Archives: June 2019
院 認識情報解析
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)
認知情報解析 実習問題
# 修正済み
# 時間がかかる場合は、繰り返し回数を減らしてください。
# 効率が悪いと思うので、適宜変更してください。
temp.state = expand.grid(loc1 = 0:2,loc2=0:2,loc3=0:2,
loc4 = 0:2,loc5=0:2,loc6=0:2,
loc7 = 0:2,loc8=0:2,loc9=0:2)
temp.state = = expand.grid(rep(list(0:2),9))
n.ones = rowSums(temp.state == 1 )
n.twos = rowSums(temp.state == 2 )
omitTwo = which(n.ones < n.twos)
omitOne = which((n.ones-1 ) > n.twos)
omitUniq = unique(c(omitOne, omitTwo))
state = temp.state[-omitUniq,]
poss.act = apply(state, 1, function(x) which(x==0))
temp.win = matrix(1:9,3)
win.idx = matrix(c(temp.win[1,],temp.win[2,],temp.win[3,],
temp.win[,1],temp.win[,2],temp.win[,3],
diag(temp.win),
diag(temp.win[3:1,])),ncol=3,byrow=T)
idx1 = c()
idx2 = c()
idxCM = c()
for (i.win in 1:nrow(win.idx)){
idx.temp = apply(state, 1, function(x) sum(x[win.idx[i.win,]]==1)==3)
idx1 = c(idx1, which(idx.temp))
idxCM.temp = apply(state, 1, function(x) sum(x[win.idx[i.win,]]==1)==2)
idxCM = c(idxCM, which(idxCM.temp))
idx.temp = apply(state, 1, function(x) sum(x[win.idx[i.win,]]==2)==3)
idx2 = c(idx2, which(idx.temp))
}
n0=apply(state,1,function(x) length((which(x==0))))
tie = which(n0==0)
Q = list()
V = list()
rew.sum = list()
rew.count = list()
policy = list()
for (i.state in 1:nrow(state)){
Q[[i.state]] = rep(0,length(poss.act[[i.state]]))
V[[i.state]] = rep(0,length(poss.act[[i.state]]))
rew.sum[[i.state]] = rep(0,length(poss.act[[i.state]]))
rew.count[[i.state]] = rep(0,length(poss.act[[i.state]]))
policy[[i.state]] = rep(1/length(poss.act[[i.state]]),length(poss.act[[i.state]]))
}
R.W = 10
R.T = 5
R.L = -10
gamma = 1
epsilon = 0.05
eta = 1
ck.result <- function(st.idx, idx1, idx2, tie){
term = F
rew = 0
result = "not terminal"
if (match(st.idx ,idx1, nomatch = 0)!=0){
rew = R.W
term = T
result = "win"
} else if (match(st.idx ,idx2, nomatch = 0)!=0){
rew = R.L
term = T
result = "lose"
} else if (match(st.idx ,tie, nomatch = 0)!=0){
rew = R.T
term = T
result = "tie"
}
return(list(rew = rew, term = term, result = result))
}
n.rep = 10000
game.hist = rep(0,n.rep)
# main loop
for (i.rep in 1:n.rep){
st.idx = 1
term = F
state.temp = rep(0,9)
state.hist1 = c()
state.hist2 = c()
repeat {
# playing game
if (length(poss.act[[st.idx]])==1){
act1 = poss.act[[st.idx]]
} else{
p.act = exp(eta*Q[[st.idx]])/sum(exp(eta*Q[[st.idx]]))
act1 = sample(poss.act[[st.idx]],1, prob = p.act)
}
state.hist1 = rbind(state.hist1,c(st.idx, act1))
state.temp[act1] = 1
st.idx = which(apply(state, 1, function(x) sum(x==state.temp) )==9)
res = ck.result(st.idx, idx1, idx2, tie)
if (res$term == T){
rew = res$rew
break
}
p.act = exp(eta*Q[[st.idx]])/sum(exp(eta*Q[[st.idx]]))
act2 = sample(poss.act[[st.idx]],1, prob = policy[[st.idx]])
state.hist2 = rbind(state.hist2,c(st.idx, act2))
state.temp[act2] = 2
st.idx = which(apply(state, 1, function(x) sum(x==state.temp) )==9)
res = ck.result(st.idx, idx1, idx2, tie)
if (res$term == T){
rew = res$rew
break
}
}
# update Q & policy
game.hist[i.rep] = res$result!="lose"
n.st = nrow(state.hist1)
#print(res$result)
if (i.rep%%100==0){print(i.rep)}
for (i.st in 1:n.st){
act.idx = which(poss.act[[state.hist1[i.st,1]]]==state.hist1[i.st,2])
rew.sum[[state.hist1[i.st,1]]][act.idx] = rew.sum[[state.hist1[i.st,1]]][act.idx] + rew
rew.count[[state.hist1[i.st,1]]][act.idx] = rew.count[[state.hist1[i.st,1]]][act.idx] + 1
Q[[state.hist1[i.st,1]]][act.idx] = rew.sum[[state.hist1[i.st,1]]][act.idx] / rew.count[[state.hist1[i.st,1]]][act.idx]
}
}
# plotting results
library(pracma)
game.hist.smooth = movavg(game.hist, 400, type="s")
plot(game.hist.smooth,type='l')
2019 データ解析基礎論A 回帰分析2
# two sample t-test
dat <- read.csv("http://www.matsuka.info/data_folder/hwsk8-17-6.csv")
boxplot(dat$ani, dat$otouto, col = c("skyblue", "orange"),ylab = "Score",
axnt ="n")
axis(1, at = c(1,2),labels = c("Ani", "Otouto"))
t.test(dat$ani, dat$otouto, var.equal=T)
dat2<-data.frame(score=c(dat$ani,dat$otouto),
order=c(rep("ani",10),rep("otouto",10)))
plot(dat2$score~as.numeric(dat2$order), pch=20, xlab="order",
ylab="score", xlim=c(0.5,2.5), cex=2, xaxt="n")
axis(1, at = c(1,2),labels = c("Ani", "Otouto"))
dat2.lm<-lm(score~order,data=dat2)
abline(dat2.lm,col='red',lwd=3)
# one sample t-test
dat.D = dat$ani - dat$otouto
boxplot(dat.D,col="skyblue",ylab="Difference")
t.test(dat.D)
plot(dat.D~rep(1,10), pch=20, xlab="", ylab="Difference", cex=3, xaxt="n")
dat.D.lm<-lm(dat.D~1)
abline(dat.D.lm,col='red',lwd=3)
# LM with 3 or more categories
dat<-read.csv("http://www.matsuka.info/data_folder/dktb312.csv")
dat2<-data.frame(result=dat$result, method=dat$method,
c1=c(rep(0,8),rep(1,8),rep(0,16)),
c2=c(rep(0,16),rep(1,8),rep(0,8)),
c3=c(rep(0,24),rep(1,8)))
dat2.lm<-lm(result~c1+c2+c3,data=dat2)
dat.lm <- lm(result~method, dat)
dat3<-data.frame(result=dat$result,
c1=c(rep(-3,8), rep(1,24)),
c2=c(rep(0,8),rep(-2,8),rep(1,16)),
c3=c(rep(0,16),rep(-1,8),rep(1,8)))
# trend analysis
dat<-read.csv("http://www.matsuka.info/data_folder/dktb321.txt")
plot(dat$result~dat$duration,data=dat[dat$method=="method.X",])
result<-dat$result[dat$method=="method.X"]
CL=c(rep(-3,5),rep(-1,5),rep(1,5),rep(3,5))
CQ=c(rep(-1,5),rep(1,5),rep(1,5),rep(-1,5))
CC=c(rep(-3,5),rep(1,5),rep(-1,5),rep(3,5))
trend.lm<-lm(result~CL+CQ+CC)
trend.lm2 <- lm(result~duration, contrasts=list(duration = "contr.poly"),data=dat.x)
summary(trend.lm2)
# trend analysis - numeric predictor
set.seed(15)
x = runif(50,0,30)
y = x^0.5 + rnorm(50,9,0.2)
dat.xy = data.frame(y=y,x=x, xsq = x^2)
plot(x,y, xlab = "Hours Studies", ylab="score",
pch=20, cex=2, col="orange")
summary(lm(y~x+xsq, data=dat.xy))
# regression diagnostic
dat<-read.csv("http://www.matsuka.info/data_folder/tdkReg01.csv")
dat.lm01<-lm(sales~price, data=as.data.frame(scale(dat)))
plot(dat.lm01,which=1)
plot(dat.lm01,which=2)
norm.vars=rnorm(300)
qqnorm(norm.vars)
qqline(norm.vars,col='red',lwd=2)
unif.vars=runif(300)
qqnorm(unif.vars)
qqline(unif.vars,col='green',lwd=2)
plot(sort(unif.vars),sort(unif.vars))
plot(sort(norm.vars),sort(unif.vars))
par(mfrow=c(2,2))
plot(dat.lm01)
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)
院：認識情報解析 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")
}
}
認知情報解析演習a Monte Carlo 01
north=c(1:3,15,1:10)
east=2:15;east[ c(3,7,11)]=c(3,7,11)
south=c(5:15,12:14)
west=c(15,1:13);west[ c(4,8,12)]=c(4,8,12)
trM=cbind(north,east,south,west)
r=-1;P=rep(0.25,4);V = rep(0,14)max.iter = 10000;
state.count=rep(0,15)
for (i.iter in 1:max.iter){
state = sample(1:14,1)
state.seq = state
while(state!=15){
action = sample(1:4,1,prob = P)
state.seq = c(state.seq,trM[state,action])
state = trM[state,action]
}
uniq.seq = unique(state.seq)
for (i.uniq in 1:(length(uniq.seq)-1)){
first.visit = which(state.seq == uniq.seq[i.uniq])[1]
V[uniq.seq[i.uniq]] = V[uniq.seq[i.uniq]] + r*(length(state.seq)-first.visit-1)
}
state.count[uniq.seq] = state.count[uniq.seq] + 1
}
V = matrix(c(0,V/state.count[1:14],0),nrow=4)
2019データ解析基礎論a 回帰分析
dat <- read.csv("http://www.matsuka.info/data_folder/hwsk8-17-6.csv")
plot(dat$otouto,dat$ani,pch=20, cex =3,
xlab = "score of younger brother",
ylab = "score of elder brother")
dat.lm <- lm(ani~otouto, data=dat)
summary(dat.lm)
abline(dat.lm, col = 'red',lwd = 2.5)
dat<-read.csv("http://www.matsuka.info/data_folder/tdkReg01.csv")
dat.reg1<-lm(sales~material,data=dat)
dat.reg2<-lm(sales~price,data=dat)
dat.reg3<-lm(sales~design,data=dat)
dat.regALL<-lm(sales~material+price+design,data=dat)
dat.regALL<-lm(sales~.,data=dat)
dat <- read.csv("http://www.matsuka.info/data_folder/hwsk8-17-6.csv")
t.test(dat$ani, dat$otouto, var.equal=T)
dat2<-data.frame(score=c(dat$ani,dat$otouto),order=c(rep("ani",10),rep("otouto",10)))
plot(dat2$score~as.numeric(dat2$order), pch=20, xlab="order",
ylab="score", xlim=c(0.5,2.5), cex=2,xaxt="n")
axis(1,c(1,2),c("ani","otouto"))
dat2.lm<-lm(score~order,data=dat2)
abline(dat2.lm,col='red',lwd=3)
dat.D = dat$ani - dat$otouto
boxplot(dat.D,col="skyblue",ylab="Difference")
t.test(dat.D)
plot(dat.D~rep(1,10),pch=20,xlab="",ylab="Difference",cex=3)
dat.D.lm<-lm(dat.D~1)
abline(dat.D.lm,col='red',lwd=3)
院：認知情報解析
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
2019 基礎実習MA01
# random number generators
x=rnorm(n=1,mean=100,sd=15)
y=runif(n=3,min=1,max=10)
N = 10000
# N = 1000
random.data = rnorm(N, mean=0, sd=1)
hist(random.data, nclass = 50, col = "navy", xlab = "Data",
probability = T, main = "Histogram of Random Data")
# density of generated data
dens = density(random.data)
lines(dens, col = "orange", lwd = 4)
# theoretical density
x = seq(-4,4,0.1)
true.norm = dnorm(x, mean = 0, sd = 1)
lines(x,true.norm, col = "green", lty = 3, lwd = 4)
legend("topleft",c("empirical", "theoretical"), lty = c(1,3),
col = c('orange','green'),lwd=4)
# random sampling
sample(1:10,3)
sample(c(“gu”,“choki”,“pa”),1)
sample(1:10)
sample(0:1, 10, replace=T)
sample(c("Head","Tail"), 10, replace=T)
sample(c("Head","Tail"), 10, replace=T, prob=c(0.9,0.1))
# flow control
for (i_loop in 1:5){print(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")
}
v1=1633;v2=355;
repeat {
r=v1%%v2
print(paste('v1 =',v1,v2 = ',v2, remainder = ',r))
v1=v2;v2=r
if (r==0){ break}
}
counter=6
repeat{
print(counter)
counter = counter + 1
if(counter>5){break}
}
counter=6
repeat{
if(counter>5){break}
print(counter)
counter+counter+1
}
認知情報解析06/05の演習問題
演習で書いたプログラムで問題はありませんでした。
以下、結果を可視化するコマンドを書き加えました：
# initalization
p.win =0.4; p.lose = 0.6; P = c(p.lose, p.win);
R = c(rep(0,100), 1); V = rep(0,101);
gamma = 1; tol = 1e-10; counter=0
cols = c("red","skyblue",'orange','black')
par(mfrow=c(1,2))
# value iteration
repeat{
counter = counter+1
delta = 0
for (i.state in 2:100) {
v <- V[i.state]
temp.v = rep(0, (i.state - 1))
for (i.bet in 1:(i.state - 1)) {
lower.B = i.state - i.bet
upper.B = min(i.state + i.bet, 101)
temp.v[i.bet] = sum(P * (R[c(lower.B, upper.B)] + gamma * V[c(lower.B, upper.B)]))
V[i.state] = max(temp.v)
}
delta <- max(abs(v-V[i.state]), delta)
}
# plotting results
if (counter==1){
plot(V[2:100], type='l', col=cols[1], lwd=2, xlab="capital", ylab="value")
} else {
if (counter<4){
lines(V[2:100], type='l', col=cols[counter], lwd=2)
}
}
if (delta < tol){break}
}
# end of value iteration
lines(V[2:100],type='l', col=cols[4], lwd=2)
legend("topleft", c("1st sweep","2nd sweep", "3rd sweep","32nd sweep") ,col=cols, lwd=2)
# identifying optimal action
policy = rep(0,101)
for (i.state in 2:100) {
temp.v = rep(0, (i.state - 1))
for (i.bet in 1:(i.state - 1)) {
lower.B = i.state - i.bet
upper.B = min(i.state + i.bet, 101)
temp.v[i.bet] = sum(P * (R[c(lower.B, upper.B)] + gamma * V[c(lower.B, upper.B)]))
policy[i.state] = which.max(round(temp.v,4))
}
}
barplot(policy,xlab="capital",ylab="Optimal action")