# 流行のモデル
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)