N = 10
M = 10000
means = rep(0, M)
for (i_rep in 1:M){
  dat <- runif(N)
  means[i_rep] <- mean(dat)
}
hist(means, xlab = "Estmatedmeans", probability = T, 
  breaks = 30, xlim = c(0,1))
x.temp = seq(0, 1, length.out = 1000)
dens <- dnorm(x.temp, mean = 0.5, sd = (1/sqrt(12*N)))
lines(x.temp, dens, col = "red", lwd = 3)
legend("topright","theoretical desiity",lwd =3, col = "red")

X = 10
Y = 1
sumXY = X + Y

summation <- function(X,Y){
  sumXY = X + Y
  return(sumXY)
}
summation(X=1, Y=10)
summation(X=5, Y=2)

CLT_example <- function(N, M){
  means = rep(0, M)
  for (i_rep in 1:M){
    dat <- runif(N)
    means[i_rep] <- mean(dat)
  }
  hist(means, xlab = "Estmated means", probability = T, 
   breaks = 30, xlim =c(0,1))
  x.temp = seq(0, 1, length.out = 1000)
  dens <- dnorm(x.temp, mean = 0.5, sd = (1/sqrt(12*N)))
  lines(x.temp, dens, col = "red", lwd = 3)
  legend("topright","theoretical desiity",lwd =3, col = "red")
}

x0 = 1e5
y0 = 10
z0 = 0
dt = 0.01
T = 10000
alpha = 1e-5
beta = 0.1
x = y = z = rep(0,T)
x[1] = x0
y[1] = y0
z[1] = z0
for(t in 1:(T-1)){
  x[t+1] = x[t] - (alpha*x[t]*y[t])*dt
  y[t+1] = y[t] + (alpha*x[t]*y[t]-beta*y[t])*dt 
  z[t+1] = z[t] + (beta*y[t])*dt
}
plot(x=1:T, x, type = "l", col = 2, lwd = 3, ylim = c(0,x0),
  xlab = "Time", ylab = "Number of people")
lines(x=1:T, y, col = 3, lwd = 3)
lines(x=1:T, z, col = 4, lwd = 3)
legend("right", c("not infected", "infected","no longer infected"), col=2:4,lwd=3)