n.subj = 5
set.seed(10)
groupA1 = rnorm(n.subj, mean = 100, sd =15)
groupA2 = rnorm(n.subj, mean = 100, sd =15)
groupID = c(rep("A1",n.subj),rep("A2",n.subj))
dat = data.frame(score = c(groupA1,groupA2),ID = groupID)
grand.mean = mean(dat$score)
mean.A1 = mean(groupA1)
mean.A2 = mean(groupA2)
plot(dat$ID, dat$score)
points(as.numeric(dat$ID), dat$score, col=c(rep('red',n.subj),rep("blue",n.subj)),pch=20)

idx = 1:(n.subj+n.subj)
plot(idx,dat$score,col=c(rep('red',n.subj),rep("blue",n.subj)),
     pch=20,ylim=c(65,145),ylab ="score",xlab="subjID")
abline(h=grand.mean,lwd=2)
abline(h=mean.A1,lty=2,col='red',lwd=2)
abline(h=mean.A2,lty=2,col='blue',lwd=2)
legend("topright",c("overall mean", "G1 mean","G2 mean"),
       lty=c(1,2,2),col=c('black','red','blue'),lwd=2)

ss.total = sum((dat$score-grand.mean)^2)
ss.group1 = n.subj*(mean.A1-grand.mean)^2
ss.group2 = n.subj*(mean.A2-grand.mean)^2
ss.group = ss.group1+ss.group2
ss.error1 = sum((groupA1-mean.A1)^2)
ss.error2 = sum((groupA2-mean.A2)^2)
ss.error = ss.error1+ss.error2
df.group = 2 - 1
df.error = (n.subj-1)+(n.subj-1)
ms.group = ss.group/df.group
ms.error = ss.error/df.error
f.value = ms.group/ms.error
p.value = 1 - pf(f.value, df.group, df.error)
summary(aov(score~ID, data=dat))

n.rep = 1e5
n.subj = 5
f.hist = rep(0,n.rep)
for (i.rep in 1:n.rep){
  groupA1 = rnorm(n.subj, mean = 100, sd =15)
  groupA2 = rnorm(n.subj, mean = 100, sd =15)
  groupID = c(rep("A1",n.subj),rep("A2",n.subj))
  dat = data.frame(score = c(groupA1,groupA2),ID = groupID)
  grand.mean = mean(dat$score)
  mean.A1 = mean(groupA1)
  mean.A2 = mean(groupA2)
  ss.total = sum((dat$score-grand.mean)^2)
  ss.group = n.subj*(mean.A1-grand.mean)^2+n.subj*(mean.A2-grand.mean)^2
  ss.error = sum((groupA1-mean.A1)^2)+sum((groupA2-mean.A2)^2)
  df.group = 2 - 1
  df.error = (n.subj-1)+(n.subj-1)
  ms.group = ss.group/df.group
  ms.error = ss.error/df.error
  f.hist[i.rep] = ms.group/ms.error
}

hist(f.hist,1000,xlim=c(0,15),probability = T)
f.temp = seq(0,15,0.05)
f.true = df(f.temp,df.group,df.error)
lines(f.temp,f.true,col='red',lty=2,lwd=2)
thres.f = qf(0.95,df.group,df.error)
abline(v = thres.f, col='green')
prop.FP = sum(f.hist > thres.f)/n.rep

# central limit theorem
n.sample = 1e6
n.data = 10
dat = matrix(runif(n.sample*n.data), nrow = n.data)
means = colMeans(dat)
sds = apply(dat,2,sd)

hist(means,100,probability = T)
x.temp  = seq(0,1, 0.01)
true.sd = sqrt(1/12)
true.dens = dnorm(x.temp,mean=0.5, sd = true.sd/sqrt(n.data))
lines(x.temp, true.dens,col='red',lwd=2)
ci.low = qnorm(0.025,mean=0.5, sd = true.sd/sqrt(n.data))
ci.high = qnorm(0.975,mean=0.5, sd = true.sd/sqrt(n.data))
abline(v=ci.low,col='blue',lwd=2)
abline(v=ci.high,col='blue',lwd=2)

# standardizing data
z = (means-0.5)/(true.sd/sqrt(n.data))
hist(z,100,probability = T)
x.tempZ = seq(-4,4,0.01)
z.dens = dnorm(x.tempZ,mean=0, sd = 1)
lines(x.tempZ, z.dens,col='red',lwd=2)
ci.lowZ = qnorm(0.025,mean=0, sd = 1)
ci.highZ = qnorm(0.975,mean=0, sd = 1)
abline(v=ci.lowZ,col='blue',lwd=2)
abline(v=ci.highZ,col='blue',lwd=2)

# t distribution
t.temp = seq(-4,4,0.01)
t.dens = dt(t.temp, df=(n.data-1))
lines(t.temp, t.dens,col='green',lwd=2)
ci.lowT = qt(0.025,df=(n.data-1))
ci.highT = qt(0.975,df=(n.data-1))
abline(v=ci.lowT,col='cyan',lwd=2)
abline(v=ci.highT,col='cyan',lwd=2)


### cell automata
nCell=201
nGen=100
res=matrix(0,nrow=nGen,ncol=nCell)
res[1,100]=1
for (i_gen in 2:nGen) {
  for (i_cell in 2:(nCell-1)) {
    res[i_gen,i_cell]=sum(res[i_gen-1,(i_cell-1):(i_cell+1)])%%2;
  } 
　#　左端
  res[i_gen,1]=(res[i_gen-1,nCell]+res[i_gen-1,1]+res[i_gen-1,2])%%2　
　#　右端
  res[i_gen,nCell]=(res[i_gen-1,(nCell-1)]+res[i_gen-1,nCell]+res[i_gen-1,1])%%2;
} 
image(res)


# general version
dec2bin<-function(num, digits=8) {
  bin=c()
  if (num==0){
    bin=0
  } else {
    while(num!=0){
      rem=num%%2
　　   num=num%/%2
      bin=c(rem,bin)
    }
  }
  if (length(bin) < digits){
    res=matrix(0,nrow=1,ncol=digits)
    res[(digits-length(bin)+1):digits]=bin
  } else {res=bin}
  return(res)
}

transFUN<-function(st,ruleID){
  output=dec2bin(ruleID,8);
  a=matrix(c(1,1,1,1,1,0,1,0,1,1,0,0,0,1,1,0,1,0,0,0,1,0,0,0),nrow=8,byrow=T)
  newSt=output[which(apply(a,1,function(x) {all(x==st)}))]
  return(newSt)
}
ECA<-function(nCell, nGen,ruleID){
  res=matrix(0,nrow=nGen,ncol=nCell)
  res[1,ceiling(nCell/2)]=1;
  for (i_gen in 2:nGen) {
    for (i_cell in 2:(nCell-1)) {
      res[i_gen,i_cell]=transFUN(res[i_gen-1,(i_cell-1):(i_cell+1)],ruleID)
     }
   res[i_gen,1]=transFUN(c(res[i_gen-1,nCell],res[i_gen-1,1],res[i_gen-1,2]),ruleID)
   res[i_gen,nCell]=transFUN(c(res[i_gen-1,(nCell-1)],res[i_gen-1,nCell],res[i_gen-1,1]),ruleID)
  }
  return(res)
}
res<-ECA(200,100,99)
image(res)