# ケーキの分配ゲーム
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)

# 中央極限定理の実験
ss.female = rnorm(6000,mean = 24, sd = 0.5)
ss.male = rnorm(6000,mean = 26, sd = 0.5)
pop = c(ss.female,ss.male)
parm.mu = mean(pop)
parm.sd = sd(pop)

n.per.one.sample = 10;
n.rep = 100000;

samp.data = matrix(sample(pop,n.per.one.sample*n.rep,replace = T),
  ncol = n.per.one.sample)
samp.mean = rowMeans(samp.data)


par(mfrow=c(1,3))
# fig 1
hist(pop,  col='skyblue',breaks = 50, probability = T, 
  main = "Population Dist.",xlab = "shoesizd")
pop.dens<-density(pop)
lines(pop.dens,col='orange',lwd=4)

# fig 2
rand.p = sample(1:n.rep,300)
dat.dens = density(samp.data[rand.p[1],])
plot(dat.dens, col='gray',xlim=c(22,28),ylim=c(0,0.5), 
  main = "Sample Dist. (N = 10)",xlab = "shoesizd")
for (i.d in 2:300){
  dat.dens = density(samp.data[rand.p[i.d],])
  lines(dat.dens, col='gray')
}

# fig 3
mean.dens<-density(samp.mean)
hist(samp.mean, col='skyblue',breaks = 50, probability = T, 
  main = 'Dist of Sample Means',xlab = "shoesizd")
lines(mean.dens,col='red',lwd=4)
x = seq(22,27,0.01)
y = dnorm(x,parm.mu, parm.sd/sqrt(n.per.one.sample))
lines(x, y ,col='darkblue',lty=3,lwd=3)
legend("topright",c("Empirical", "Theoretical"),
  lty=c(1,3),col=c("red","blue"),lwd=4,cex=1)