options(digits=3)

# Generate 10 random normal observations with mean 80, s.d. 15
rnorm(10, 80, 15)
# Re-issue this command multiple times to see new sets of random #s

# For each set of random #s, print the mean and s.d.
x <- rnorm(10, 80, 15)
c(mean(x), sqrt(var(x)))
# Repeat several times

# Make a histogram and ecdf from a sample of size 500
x <- rnorm(500, 80, 15)
hist(x, nclass=30)
ecdf(x, q=.5)   # show median using reference lines

# Use two different program to generate "reps" samples of size 10,
# saving the sample mean from each sample
sim <- function(reps, n, mu, sd) {
   means <- single(reps)  # set aside "reps" numbers
   for(i in 1:reps) means[i] <- mean(rnorm(n, mu, sd))
   means
}

sim(5, 10, 80, 15)   # repeat several times

# The second method uses sapply, which ordinarily operates on all
# variables in a list.  Here there's no list but we want to do a
# separate calculation for each of the "reps" elements

sapply(1:5, function(...) mean(rnorm(10, 80, 15)))

# Use this method to draw a sample of 500 means, each from a
# sample of size 10 with mean 80 sd 15

means <- sapply(1:500, function(...) mean(rnorm(10, 80, 15)))
# Print first 5 of these means
means[1:5]

# Compute mean, var, s.d. of the 500 means
mean(means); var(means); sqrt(var(means))

# Compare variance with (15^2) / 10, s.d. with 15/sqrt(10)
# 10 = individual sample sizes
(15^2) / 10; 15/sqrt(10)

# Show distribution of these 500 sample means
hist(means, nclass=30); ecdf(means, q=.5)


#--------------------------------------------------------------------

# Compute probability of getting 0, 1, 2, 3, 4 heads out of 4 tosses
dbinom(0:4, 4, .5)
plot(0:4, dbinom(0:4, 4, .5))

# Two methods to get Prob[<= x heads], x=0-4
cumsum(dbinom(0:4, 4, .5))
pbinom(0:4, 4, .5)

# Get expected # head out of 4 tosses
sum(0:4 * dbinom(0:4, 4, .5))

# Plot probability of getting 0:20 heads out of 20 tosses
plot(0:20, dbinom(0:20, 20, .5))

# Flip a coin 1000 times and observe # of heads
rbinom(1, 1000, .5)

# Do this 2000 times and plot # heads vs. toss number
plot(1:2000, heads <- rbinom(2000, 1000, .5))

# Compute mean and variance of the 2000 realizations of the
# binomial random variable from n=1000 p=.5; compare to theoreticals
mean(heads); var(heads); 1000*.5; 1000*.5*(1-.5)