# Code for example done in class on January 17, 2019

# Code to examine the Central Limit Theorem for 
# a single sample size

n <- 20
N <- 1000
Z <- rep(NA, N)
for (i in 1:N) {
  # Generate a sample from an exponential distribution.
  x <- rexp(n, 1) 
  # Calculate the standardized mean
  Z[i] <- sqrt(n)*(mean(x)-1)
}
hist (Z, breaks = 40, prob=T)
lines(seq(-4,4,by=0.01), dnorm(seq(-4,4,by=0.01)),
      col="red")

# An alternative using implicit loops

# Generate all of the data at once and store in a matrix
x.data <- matrix(rexp(n*N,1), nrow=N)
Z <- sqrt(n)*(rowMeans(x.data)-1)
hist (Z, breaks = 40, prob=T)
lines(seq(-4,4,by=0.01), dnorm(seq(-4,4,by=0.01)),
      col="red")
# Aternatively we could use 
# apply(x.data, 1, mean) in place of rowMeans(x.data)

# Extend this to a number of sample sizes.
# We use two nested loops here
nn <- c(3, 10, 20, 50, 100)
N <- 1000
output <- matrix(NA, nrow=N, ncol=length(nn))
j <- 1
for (n in nn) {
  Z <- rep(NA, N)
  for (i in 1:N) {
    # Generate a sample from an exponential distribution.
    x <- rexp(n, 1) 
    # Calculate the standardized mean
    Z[i] <- sqrt(n)*(mean(x)-1)
  }
  # Store the output in the jth column and increase j
  output[,j] <- Z
  j <- j+1
}


# Examle of a While loop
# Newton-Raphson Algorithm to solve f(x)=0
f <- function(x)
  # The function for which we wish to find the root
  3*x^4-2*x^3-3*x^2-2*x+3

# Plot the function and ensure there is root.
# Also get its approximate value as an initial point.
curve(f)
abline(h=0)

f.prime <- function(x) 
  # The derivative of the function f
  12*x^3-6*x^2-6*x-2

x <- 0.7
continue <- T
R <- 100
eps <- 1e-12
output <- matrix(NA, nrow=R+1, ncol=3)
i <- 1
while (continue) {
   output[i,] <- c(x, f(x), f.prime(x)) 
   x <- x-output[i,2]/output[i,3]
   i <- i+1
   continue <- (i <= R)&(abs(f(x)) >= eps)
}
# Remove rows of output not used  
output <- output[!is.na(output[,1]),]
solution <- output[nrow(output),]