# Simulate random variables in R and plot a histogram, for Assignment 3



# Example 1, Make a list of standard normal random variables

n = 5           # number of random numbers
x = rnorm(n)     # get a list of that many standard normals

cat("\nExample 1, a list of ", n, " random numbers\n\n")
for (k in 1:n){
   cat("     x[", k , "] = ", x[k],"\n")
  }

cat("\n x is: ", x, "\n")



#  Example 2, make a list the slow way

cat("\nExample 2, another list of ", n, " random numbers\n\n")
n = 5
x = rnorm(1)
for ( k in 2:n){
   x_k = rnorm(1)[1]
   x   = c( x, x_k)
   cat("   after appending ", x_k, ", x is ", x, "\n")
  }



# Example 3, Make n more, to see the seed at work

n = 5           # number of random numbers
s = 1
set.seed(s)
x = rnorm(n)     # get a list of that many standard normals

cat("\nExample 3, a list of ", n, " random numbers with seed = ", s, "\n\n")
cat("   The list x is ", x, "\n")



# Example 4, Make make n more, again, to see the seed at work

n = 5           # number of random numbers
s = 1
set.seed(s)
x = rnorm(n)     # get a list of that many standard normals

cat("\nExample 4, a list of ", n, " random numbers with seed = ", s, "\n\n")
cat("   The list x is ", x, "\n")



# Example 5, Make n more, yet again, to see the seed at work

n = 5           # number of random numbers
s = 2
set.seed(s)
x = rnorm(n)     # get a list of that many standard normals

cat("\nExample 5, a list of ", n, " random numbers with seed = ", s, "\n\n")
cat("   The list x is ", x, "\n")





# Example 6, make a histogram

n = 1000000
x = rnorm(n)
Title = sprintf("%7d standard normals, the popup", n)

#    Make a popup histogram

hist( x,
      breaks      = 50,          #  The number of histogram bins
      main        = Title,       #  Add the title
      probability = TRUE)        #  Plot the estimated probability density

#   Put the same histogram into a pdf file

pdf("NormalHistogram.pdf")
Title = sprintf("%7d standard normals, the .pdf file", n)
hist( x,
      breaks      = 50,          #  The number of histogram bins
      main        = Title,       #  Add the title
      probability = TRUE)        #  Plot the estimated probability density
dev.off()

f_0 = 1/sqrt(2*pi)             # this is 1/sqrt(2*pi)e^(-x^2/2) with x = 0.
cat("\nExample 6, the probability density at x = 0 should be ", f_0,"\n")






# Example 7, make a histogram of the max of normals

n = 100000         # number of experiments
L = 2              # number of normals to take the max of

x = c(1:n)          # pre-allocate the results vector x, for efficiency

for (k in 1:n){
   y    = rnorm(L)  # Generate the L normals
   x[k] = max(y)    # Find the max, record the result of experiment k
  }

Title = sprintf("Maximum of %5d normals, %7d samples", L, n)

#    Make a popup histogram

dev.new()                        #  Open this plot in a different popup so you can see both
hist( x,
      breaks      = 50,          #  The number of histogram bins
      main        = Title,       #  Add the title
      probability = TRUE)        #  Plot the estimated probability density

#   Put the same histogram into a pdf file

pdf("ExperimentHistogram.pdf")
hist( x,
      breaks      = 50,          #  The number of histogram bins
      main        = Title,       #  Add the title
      probability = TRUE)        #  Plot the estimated probability density
dev.off()



#  Example 8, exponential random variables.

n   = 1000
lam = .5

T = 1:n          # Create a list without caring what's in it
U = runif(n)
for (k in 1:n){
   T[k] = -(1/lam) * log( U[k])
  }

N_1 = 0      # Count the number of T_k > t_1
N_2 = 0
t_1 = 1
t_2 = 3

for (k in 1:n){
   if ( T[k] > t_1) {
      N_1 = N_1 + 1
     }
   if ( T[k] > t_2) {
      N_2 = N_2 + 1
     }
  }

P1_em = N_1/n
P2_em = N_2/n
P1_th = exp( -lam*t_1)
P2_th = exp( -lam*t_2)

cat("\nExample 8: exponential random variables\n\n")

cat("    Probability of T >", t_1, ",   empirical = ", P1_em, ", and theoretical = ", P1_th, "\n")
cat("    Probability of T >", t_2, ",   empirical = ", P2_em, ", and theoretical = ", P2_th, "\n")

Title = sprintf("Exponential with rate = %5.3f, %7d samples", lam, n)

dev.new()
hist( T,
      breaks      = 50,          #  The number of histogram bins
      main        = Title,       #  Add the title
      probability = TRUE)        #  Plot the estimated probability density

pdf("ExponentialHistogram.pdf")
hist( T,
      breaks      = 50,          #  The number of histogram bins
      main        = Title,       #  Add the title
      probability = TRUE)        #  Plot the estimated probability density
dev.off()