認知情報解析　01

# 流行のモデル
ryuko <- function(x0, y0, z0, alpha, beta, T){
dt = 0.01
total.p = x0+y0+z0
x = rep(0,T);y = rep(0,T);z = rep(0,T)
x[1]=x0;y[1]=y0;z[1]=z0;
for (i.t in 1:(T-1)){
x[i.t + 1] =x[i.t]-(alpha*x[i.t]*y[i.t])*dt
y[i.t + 1] =y[i.t]+(alpha*x[i.t]*y[i.t]-beta*y[i.t])*dt
z[i.t + 1] =z[i.t]+(beta*y[i.t])*dt
}
plot(1:T,x/total.p,type='l',col='red',xlab = "Time",
ylab = "Proportion", lwd = 2)
lines(1:T, y/total.p, col='blue', lwd = 2)
lines(1:T, z/total.p, col='green', lwd = 2)
return(data.frame(x=x,y=y,z=z))
}
# 大数の法則に関するシミュレーション

# サイコロをN回ふって６が出る確率の推移を検証する関数
die.6 <- function(N){
die <- sample(1:6, N, replace=T)
six.cumsum = cumsum(die==6)
six.prob = six.cumsum/(1:N)
return(six.prob=six.prob)
}

# 関数の実行＆結果の可視化
N = 3000
result = die.6(N)
plot(1:N, result, ylim=c(0,1), col='red', lwd=3, type='l',
xlab="Trial",ylab="Prob.")
abline(h = 1/6, lty=2, lwd=2)
legend("topright",c("Theoretical","Empirical"),lty=c(2,1),
col=c("black","red"), lwd=2)

# 大数の法則を繰り返し検証するスクリプト
N = 3000; n.rep = 500
result = matrix(0, nrow=n.rep, ncol=N)
plot(1,1,type='n',xlim = c(0,N), ylim=c(0,1),
xlab="Trial", ylab="P(die = 6)")
for (i.rep in 1:n.rep){
result[i.rep,] = die(N)
lines(result[i.rep,],col='gray')
}
lines(1:N, colMeans(six.p), col='red', lwd = 2)
abline(h = 1/6, lty=2, lwd=2)
legend("topright",c("Theoretical","Average","Indiv. trial"),lty=c(2,1,1),
col=c("black","red","gray"), lwd=2)

# 大数の法則を繰り返し検証するスクリプトを関数にしたもの
die.6Mult <- function(N, n.rep) {
plot(1,1,type='n',xlim = c(0,N), ylim=c(0,1), xlab="Trial",
ylab="P(die = 6)")
die <- matrix(sample(1:6, N*n.rep, replace=T), nrow = n.rep)
six.cumsum = apply(die,1, function(x){cumsum(x==6)})
six.prob = six.cumsum/(1:N)
for (i.rep in 1:n.rep){
lines(six.prob[,i.rep],col='gray')
}
lines(1:N, rowMeans(six.prob), col='red', lwd = 2)
abline(h = 1/6, lty=2, lwd=2)
legend("topright",c("Theoretical","Average","Indiv. trial"),lty=c(2,1,1),
col=c("black","red","gray"), lwd=2)
return(six.prob=six.prob)
}
#  上記の関数の実行例
result <- die.6Mult(3000,500)

2018 DAA W03 Data Viz & Descrptive statistics

par(mfrow=c(1,2))
hist(dat[dat\$gender=='F',]\$h, main="Dist. of Height for Female Participants",
xlab="Height", xlim=c(140,190), probability=T)
dens.F = density(dat[dat\$gender=='F',]\$h)
lines(dens.F, col='blue',lwd=2)
hist(dat[dat\$gender=="M",]\$h,
main="Dist. of Height for Male Participants",
xlab="Height", xlim=c(140,190), probability=T,ylim=c(0,0.08))
dens.M = density(dat[dat\$gender=='M',]\$h)
lines(dens.M, col='green', lwd=2)

par(mfrow=c(1,1))
plot(dens.F,col='blue',lwd=2, ylab='density', xlim=c(140,190),
main="Dist. of Height by gender",xlab='Height')
lines(dens.M,col='green',lwd=2)
legend("topleft", c('Female','Male'), col=c('blue','green'), cex=1.5,lwd=2)

text(x = 157.5, y = 0.04, 'Female', col='blue', cex=2)
text(x = 170, y = 0.04,'Male', col='green', cex=2)

plot(dat\$shoesize, dat\$h, main="Relationship b/w shoesize and height",
xlab = "shoesize", ylab="height", pch=19, col="red")

txt = paste("r =",round(cor(dat\$shoesize,dat\$h), 4))
text(22, 175, txt, cex = 1.5)

abline(h = mean(dat\$h), col='blue');
abline(v = mean(dat\$shoesize), col='green');
abline(lm(dat\$h~dat\$shoesize), lty=2, lwd=2)

plot(dat[dat\$gender=='F',]\$shoesize, dat[dat\$gender=='F',]\$h,
main="Relationship b/w shoesize and height", xlab='shoesize', ylab='height',
cex.lab=1.5, pch=19, col='blue', xlim=c(20,29), ylim=c(140,190))
points(dat[dat\$gender=='M',]\$shoesize,dat[dat\$gender=='M',]\$h,
pch = 15, col = 'green')
legend("topleft", c('Female','Male'), pch =c(19,15),
col = c('blue','green'), cex = 1.5)

plot(dat, pch=20, col='blue')

plot(dat.pca, pch = rownames(dat.pca), cex = 1.7, col = 'blue')

2018 DAA W02 Data visualization

v1 = seq(-3,3,0.1)
v2 = v1^2
plot(x = v1, y = v2)
plot(v1, v2, col = 'red')
plot(v1, v2, col=“red”, pch = 20)
plot(v1, v2, col="red", pch = 20, cex = 3)
plot(v1, v2, type = "l")
plot(v1, v2,type="l",lty=4,lwd=3)

plot(v1, v2, main = "THIS IS THE TITLE",
xlab = "Label for X-axis",
ylab = "Label for Y-axis")

plot(v1, v2, main = "THIS IS THE TITLE", cex.lab = 1.5,
xlab = "Label for X-axis",ylab = "Label for Y-axis")

plot(v1, v2, main = "TITLE", xlab = "X here",ylab = "Y here",
xlim = c(-3.5, 3.5),
ylim = c(-0.5, 10))

plot(v1, v2, col = "blue", type = "o", lty = 2, pch = 19,
cex.lab = 1.5, lwd = 3, main = "Y=X*X", xlab = "X",
ylab="X*X", xlim=c(-3.5,3.5), ylim=c(-0.5, 10))

hist(dat\$h)
hist(dat\$h, breaks = 20, main = "Histogram of Height",
xlab = "Height", col = 'blue', xlim = c(140, 190))

dens<-density(dat\$h); # 確率密度の算出
hist(dat\$h, main = "Histogram of Height", xlab = "Height",
xlim = c(140,190), probability = T)
lines(dens, lwd = 2, col = "red", lty=2)

plot(v1, v2, col = "blue", type = "l", pch = 19, cex.lab = 1.5, lwd = 3, xlab = "X",
ylab="f(X)", xlim=c(-3.5,3.5), ylim=c(-0.5, 10))
lines(v1, v1^3, col='red',lwd = 3)
legend("bottomright", c("x^2","x^3"), col=c('blue','red'), lwd=2)

boxplot(dat\$h,main="Boxplot of Height", ylab="Height", col='cyan', ylim=c(140,190))
boxplot(dat\$h,main="Boxplot of Height", xlab="Height", col='orange', horizontal=T)

boxplot(dat\$h ~ dat\$gender,
main="Distribution of Height by Gender",
ylab="Gender", xlab="Height", col=c('blue','cyan'),
ylim=c(140,190), horizontal=T)

boxplot(dat\$h ~ dat\$gender + dat\$affil,
main="Distribution of  Height by Gender and Affiliation",
ylab="Gender x Affiliation", xlab="Height", col=c("blue","cyan","red","magenta"),
ylim=c(140,190),horizontal=T)

interaction.plot(dat\$gender,
dat\$affil,
dat\$h,
pch=c(20,20),
col=c("skyblue","orange"),
xlab="gender", ylab="height",
lwd=3,type='b',cex=2,
trace.label="Affiliation")

par(mfrow=c(1,2))
hist(dat[dat\$gender=="F",]\$h, main="Dist. of Height for Female Participants", xlab="Height", xlim=c(140,190), probability=T)
dens.F = density(dat[dat\$gender=='F',]\$h)
lines(dens.F, col='blue',lwd=2)
hist(dat[dat\$gender=="M",]\$h, main="Dist. of Height for Male
Participants", xlab="Height", xlim=c(140,190), probability=T,ylim=c(0,0.08))
dens.M = density(dat[dat\$gender=='M',]\$h)
lines(dens.M, col='green', lwd=2)

par(mfrow=c(1,1))
plot(dens.F,col='blue',lwd=2, ylab='density', xlim=c(140,190), main="Dist. of Height by gender",xlab='Height')
lines(dens.M,col='green',lwd=2)
legend("topleft", c('Female','Male'), col=c('blue','green'), cex=1.2,lwd=2)

plot(dat\$shoesize, dat\$h, main="Relationship b/w shoesize and height",
xlab = "shoesize", ylab="height", pch=19, col="red")