Monthly Archives: October 2017

Posted on by

ケーキの分配ゲーム

説明はここにあります。

# ケーキの総数
n.cake = 10 

# 利得行列
pay.mat = matrix(0,(n.cake+1),(n.cake+1))
for (i.cake in 1:n.cake){
  pay.mat[(i.cake+1),1:(n.cake-i.cake+1)] =i.cake
}

# 初期化（各戦略の確率や時間など）
p.cake = runif(n.cake+1)
p.cake = p.cake/sum(p.cake)
max.time = 50
dt = 0.01
t = seq(0,max.time,dt)
n.iter = length(t)
p.hist = matrix(0,nrow = n.iter, ncol = (n.cake+1))
p.hist[1,] = p.cake

# メインのスクリプト
for (i.time in 2:n.iter){
  W = colSums(p.cake*t(pay.mat))
  W.ave = sum(outer(p.cake,p.cake)*pay.mat)
  p.cake = p.cake + p.cake*(W - W.ave)/W.ave * dt
  p.hist[i.time,] = p.cake
}

# 結果の可視化
plot(p.hist[,1],ylim=c(0,1),type='l',lwd=2,ylab = 'Proportion',xlab="time")  
for (i.strat in 2:(n.cake+1)){
  lines(p.hist[,i.strat],col=i.strat,lwd=2)
}
legend("topleft",paste("request = ",0:10),col=1:(n.cake+1),lwd =2,cex=0.75)
Posted on by

データ解析基礎論B W03

# RB
dat<-read.csv("http://www.matsuka.info/data_folder/dktb316.txt")
colnames(dat)<-c("id",'method','subj','words')
dat.aov<-aov(words~method+subj+Error(subj:method),data=dat)
dat.aov.sum = summary(dat.aov)

# Tukey HSD
mse <- dat.aov.sum[[1]][[1]][3,3]
q<-qtukey(0.95,4,df=21)
hsd<-q*sqrt(mse/8)
dat.means=tapply(dat$words, dat$method, mean)
ck.hsd<-abs(outer(dat.means,dat.means,'-'))>hsd

# SPF
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")

source("http://peach.l.chiba-u.ac.jp/course_folder/tsme.txt")
dat<-read.csv("http://www.matsuka.info/data_folder/dktb3211.txt")
dat.aov<-aov(result~method*duration+Error(s+s:duration),dat)
Posted on by

データ解析基礎論B　ー　期末テストの解説

# Q2 example
# visualization
interaction.plot(dat.Q2$med,dat.Q2$ingred,dat.Q2$bp,type='o',
   col = c("skyblue", "coral"),lwd =3, cex = 2,
   pch=c(18,19),xlab = "quantitily", ylab = "Blood Pressure",
   trace.label="type of medicine",legend = T)
# descriptive stats
tapply(dat.Q2$bp,list(dat.Q2$med,dat.Q2$ingred),mean)
# ANOVA
dat.Q2.aov = aov(bp~med*ingred,data=dat.Q2)
summary(dat.Q2.aov)
# Linear Regression
dat.Q2.lm = lm(bp~med*ingred, contrasts=list(med=contr.poly), data=dat.Q2)
summary(dat.Q2.lm)
anova(dat.Q2.lm)

# Q3 example 1
# visualization
plot(hours~btype, data= dat.Q3V1)
# ANOVA
dat.Q3V1.aov = aov(hours~btype,data =dat.Q3V1)
summary(dat.Q3V1.aov)
TukeyHSD(dat.Q3V1.aov)
# Linear Regression
dat.Q3V1.lm = lm(hours~btype, data = dat.Q3V1)
summary(dat.Q3V1.aov)

# Q3 example 2
# visualization
plot(dat.Q3V1$hours, dat.Q3V1$A,type = 'p',pch=20,cex =2, 
   xlab = "hours cleaned", ylab= "Blood Type",yaxt ="n")
axis(2, at = c(0,1), labels = c("not A", "A"))
# GLM
dat.Q3V1.glm=glm(A~hours,family="binomial",data=dat.Q3V1 )
# plotting result
x.temp = seq(100,330,1)
glm.cf = coef(dat.Q3V1.glm)
y.temp = 1/(1+exp(-1*(glm.cf[1]+glm.cf[2]*x.temp)))
lines(x.temp,y.temp,col='red',lwd=2)

# Q4 example
plot(score~group,data=dat.Q4)
dat.Q4.G = lm(score~group,data=dat.Q4)

#
dat.Q4.NG1 = lm(score~income, data=dat.Q4)
dat.Q4.NG2 = lm(score~gender, data=dat.Q4)
dat.Q4.NG3 = lm(score~study, data=dat.Q4)
dat.Q4.NG4 = lm(score~height, data=dat.Q4)

dat.Q4.NG5 = lm(score~study+group,data=dat.Q4)
dat.Q4.NG6 = lm(score~study+income,data=dat.Q4)