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)