7  Simulations

7.1  The central limit theorem

> n=10;p=.25;S= rbinom(1,n,p)
> (S - n*p)/sqrt(n*p*(1-p))
[1] -0.3651484
  

> n = 10;p = .25;S = rbinom(100,n,p)
> X = (S - n*p)/sqrt(n*p*(1-p)) # has 100 random numbers
  

> hist(X,prob=T)    
  

---

Figure 31: Scaled binomial data is approximately normal(0,1)

---

7.2  For loops

> results =numeric(0)          # a place to store the results
> for (i in 1:100) {           # the for loop
+ S = rbinom(1,n,p)            # just 1 this time
+ results[i]=(S- n*p)/sqrt(n*p*(1-p)) # store the answer
+ }

> primes=c(2,3,5,7,11);
## loop over indices of primes with this
> for(i in 1:5) print(primes[i])
## or better, loop directly
> for(i in primes) print(i)
  

> results = c();
> mu = 0; sigma = 1
> for(i in 1:200) {
+ X = rnorm(100,mu,sigma)       # generate random data
+ results[i] = (mean(X) - mu)/(sigma/sqrt(100))
+ }
> hist(results,prob=T)

7.3  Normal plots

> x = rnorm(100,0,1);qqnorm(x,main='normal(0,1)');qqline(x)
> x = rnorm(100,10,15);qqnorm(x,main='normal(10,15)');qqline(x)
> x = rexp(100,1/10);qqnorm(x,main='exponential mu=10');qqline(x)
> x = runif(100,0,1);qqnorm(x,main='unif(0,1)');qqline(x)
  

---

Figure 33: Some normal plots

---

7.4  Using simple.sim and functions

> f = function () {
+ S = rbinom(1,n,p)           
+ (S- n*p)/sqrt(n*p*(1-p))
+ }

> x=simple.sim(100,f)
> hist(x)

> f = function(n=100,p=.5) {
+ S = rbinom(1,n,p)           
+ (S- n*p)/sqrt(n*p*(1-p))
+ }

> simple.sim(1000,f,100,.5)

• Functions have a special keyword function as in

  
  > the.range = function (x) max(x) - min(x)    
    

  which returns the range of the vector x. (Already available with range.) This tells R that the.range is a function, and its arguments are in the braces. In this case (x).
• If a function is a little more complicated and requires multiple commands you use braces (like a for loop). The last value computed is returned. This example finds the IQR based on the lower and upper hinges and not the quantiles. It uses the results of the fivenum command to get the hinges

  
  > find.IQR = function(x) {
  + five.num = fivenum(x)         # for Tukey's summary
  + five.num[4] - five.num[2]
  + }    
    

  The plus sign indicates a new line and is generated by R -- you do not need to type it. (The five number summary is 5 numbers: the minimum, the lower hinges, the median, the upper hinge, and the maximum. This function subtracts the second from the fourth.)
• A function is called by its name and with parentheses. For example

  
  > x = rnorm(100)               # some sample data
  > find.IQR                     # oops! no argument. Prints definition.
  function(x) {
  five.num = fivenum(x)
  five.num[4] - five.num[2]
  }
  > find.IQR(x)                  # this is better
  [1] 1.539286
    

> f = function(n=100,mu=0,sigma=1) {
+ nos = rnorm(n,mu,sigma) 
+ (mean(nos)-mu)/(sigma/sqrt(n))
+ }
  

> simulations = simple.sim(100,f,100,5,5) 
> hist(simulations,breaks=10,prob=TRUE)

  > f = function(n=100,mu=10) (mean(rexp(n,1/mu))-mu)/(mu/sqrt(n))

> xvals = seq(-3,3,.01) # for the density plot 
> hist(simple.sim(100,f,1,10),probability=TRUE,main="n=1",col=gray(.95))
> points(xvals,dnorm(xvals,0,1),type="l") # plot normal curve 
  ... repeat for n=5,15,50

7.5  Problems

7.1

Do two simulations of a binomial number with n=100 and p=.5. Do you get the same results each time? What is different? What is similar?

7.2

Do a simulation of the normal two times. Once with n=10, µ=10 and s=10, the other with n=10, µ = 100 and s=100. How are they different? How are they similar? Are both approximately normal?

7.3

The Bernoulli example is also skewed when p is not .5. Do an example with n=100 and p=.25, p=.05 and p=.01. Is the data approximately normal in each case? The rule of thumb is that it will be approximately normal when np³5 and n(1-p)³5. Does this hold?

7.4

The normal plot is a fancy way of checking if the distribution looks normal. A more primitive one is to check the rule of thumb that 68% of the data is 1 standard deviation from the mean, 95% within 2 standard deviations and 99.8% within 3 standard deviations.

Create 100 random numbers when the X_i are normal with mean 0 and standard deviation 1. What percent are within 1 standard deviation of the the mean? Two standard deviations, 3 standard deviations? Is your data consistent with the normal?

(Hint: The data is supposed to have mean 0 and variance 1. To check for 1 standard deviation we can do

> k = 1;sigma = 1
> n = length(x)
> sum( -k*sigma <x & x< k*sigma)/n

Read the & as "and" and this reads as -- after simplification-- ``-1 less than x and x less than 1''. This is the same as P(-1<x<1).)

7.5

It is interesting to graph the distribution of the standardized average as n increases. Do this when the X_i are uniform on [0,1]. Look at the histogram when n is 1, 5, 10 and 25. Do you see the normal curve taking shape?

(A rule of thumb is that if the X_i are not too skewed, then n>25 should make the average approximately normal. You might want

 
> f=function(n,a=0,b=1) {
    mu=(b+a)/2
    sigma=(b-a)/sqrt(12)
    (mean(runif(n,a,b))-mu)/(sigma/sqrt(n))
}

where the formulas for the mean and standard deviation are given. )

7.6

A home movie can be made by automatically showing a sequence of graphs. The system function System.sleep can insert a pause between frames. This will show a histogram of the sampling distribution for increasingly large n

> for (n in 1:50) {
+ results = c()
+ mu = 10;sigma = mu
+ for(i in 1:200) {
+ X = rexp(200,1/mu)
+ results[i] = (mean(X)-mu)/(sigma/sqrt(n))
+ }
+ hist(results)
+ Sys.sleep(.1)
+ }
  

Run this code and take a look at the movie. To rerun, you can save these lines into a function or simply use the up arrow to recall the previous set of lines. What do you see?

7.7

Make normal graphs for the following random distributions. Which of them (if any) are approximately normal?

1. rt(100,4)
2. rt(100,50)
3. rchisq(100,4)
4. rchisq(100,50)

7.8

The bootstrap technique simulates based on sampling from the data. For example, the following function will find the median of a bootstrap sample.

> bootstrap=function(data,n=length(data)) {
+  boot.sample=sample(data,n,replace=TRUE)
+  median(boot.sample)
+ }

Let the data be from the built in data set faithful. What does the distribution of the bootstrap for the median look like? Is it normal? Use the command:

> simple.sim(100,bootstrap,faithful[['eruptions']])  

7.9

Depending on the type of data, there are advantages to the mean or the median. Here is one way to compare the two when the data is normally distributed

> res.median=c();res.mean=c() # initialize
> for(i in 1:200) {             # create 200 random samples
+ X = rnorm(200,0,1)
+ res.median[i] = median(X);res.mean[i] = mean(X)
+ }
> boxplot(res.mean,res.median)  # compare
  

Run this code. What are the differences? Try, the same experiment with a long tailed distribution such as X = rt(200,2). Is there a difference? Explain.

7.10

In mathematical statistics, there are many possible estimates for the center of a data set. To choose between them, the one with the smallest variance is often taken. This variance depends upon the population distribution. Here we investigate the ratio of the variances for the mean and the median for different distributions. For normal(0,1) data we can check with

> median.normal = function(n=100) median(rnorm(100,0,1))    
> mean.normal = function(n=100) mean(rnorm(100,0,1))
> var(simple.sim(100,mean.normal)) /
+ var(simple.sim(100,median.normal))
[1] 0.8630587
  

The answer is a random number which will usually be less than 1. This says that usually the variance of the mean is less than the variance of the median for normal data. Repeat using the exponential instead of the normal. For example:

> mean.exp  =  function(n=100) mean(rexp(n,1/10))
> median.exp = function(n=100) median(rexp(n,1/10))
  

and the t-distribution with 2 degrees of freedom

> mean.t = function(n=100) mean(rt(n,2))
> median.t = function(n=100) median(rt(n,2))
  

Is the mean always better than the median? You may also find that side-by-side boxplots of the results of simple.sim are informative.