Simulation studies using R

Common continuous distributions include

Exponential

(exp) The parameter is rate (defaults to 1). The mean of the distribution is 1/rate.

Normal (or Gaussian)

(norm) The most famous distribution in statistics, this is the well-known "bell-curve". Parameters of the distribution are mean (defaults to 0) and sd (defaults to 1).

Uniform

(unif) Parameters are min (defaults to 0) and max (defaults to 1).

Common discrete distributions include

Binomial

(binom) Parameters are size, the number of trials, and prob, the probability of success on a single trial (no defaults).

Geometric

(geom) Parameter is prob, the probability of success on each independent trial. Note that the distribution is defined in terms of the number of failures before the first success.

Poisson

(pois) The parameter is lambda, the mean.

Less common continuous distributions include

Beta

(beta) Parameters are shape1 and shape2 (without defaults) and ncp, the non-centrality parameter, with a default of 0, corresponding to the central beta.

Cauchy

(cauchy) Parameters are location (defaults to 0) and scale (defaults to 1).

Gamma

(gamma) Parameters are shape and one of rate or scale. The last two parameters are inverses of one another and both default to 1.

Logistic

(logis) Parameters are location and scale.

Log-normal

(lnorm) Parameters are meanlog and sdlog.

Weibull

(weibull) Parameters are shape and scale.

Continuous distributions describing sample statistics from a Gaussian population include

Chi-square

(chisq) Parameters are df and ncp, the non-centrality parameter.

F

(f) Parameters are df1 and df2, the numerator and denominator degrees of freedom. An optional parameter is ncp, the non-centrality parameter. If omitted a central F distribution is assumed.

t

(t) Parameters are df and ncp, the non-centrality parameter. If omitted a central t distribution is assumed.

Less common discete distributions include the hypergeometric (hyper) and the negative binomial (nbinom).

The random number stream depends on a seed value. If we want to produce a reproducible example (so, for example, we can talk about properties of a particular sample and the reader can generate the same sample for herself so she can examine it) then we set the seed to a known value before generating the sample. The stored seed is a complicated structure but we can set it to an integer, often something trivial like

set.seed(1)                             # Ensure a reproducible stream

If we simulate a sample,

(s1 <- rnorm(5))

[1] -0.6264538  0.1836433 -0.8356286  1.5952808  0.3295078

then similate a second sample

(s2 <- rnorm(5))

[1] -0.8204684  0.4874291  0.7383247  0.5757814 -0.3053884

we get different values. However, if we reset the seed to 1 and then simulate the second sample we reproduce the original values.

set.seed(1)
(s2 <- rnorm(5))

[1] -0.6264538  0.1836433 -0.8356286  1.5952808  0.3295078

We can see from the output that the printed values of the two samples are similar. However, the printed values are rounded. A more reliable check is

all.equal(s1, s2)

[1] TRUE

The theoretical expected value and variance of a distribution are defined in terms of the parameters. Sometimes the definitions are simple and easy to remember; other times they aren't. A simple way to examine how, say, the expected value depends on the parameter is to simulate a very large sample from the distribution and evaluate the sample mean.

Recall that the parameter rate in the R d-p-q-r functions for the exponential distribution is the inverse of the expected value. We would anticipate that the sample mean for a large sample will be very close to the theoretical mean, 1/rate.

mean(rexp(1e6, rate = 1))

[1] 1.000786

mean(rexp(1e6, rate = 0.5))

[1] 2.000935

mean(rexp(1e6, rate = 2))

[1] 0.5007178

So the theoretical value of 1/rate seems to hold.

What about the variance or the standard deviation. I can remember that one of them is also 1/rate but sometimes I forget which one.

var(rexp(1e6, rate = 1))

[1] 1.000989

var(rexp(1e6, rate = 0.5))

[1] 4.002376

var(rexp(1e6, rate = 2))

[1] 0.2494673

So it looks as if the variance is 1/rate² which would mean that it is the standard deviation that is 1/rate.

• Exercises
  
  1. Evaluate the sample mean, variance and standard deviation of one or more large samples from the Geometric distribution for different values of prob (remember that this parameter represents a probability and must be between 0 and 1). Look at the Wikipedia description for the Geometric distribution, which describes two ways of writing the distribution
     • the number of failures before the first success
     • the number of trials until the first success

     Which version is implemented in R?

  2. Evaluate the sample mean, variance and standard deviation of one or more large samples from the Uniform distribution for different values of the parameters min and max. How does the variance related to min and max.
  3. Evaluate the sample mean, variance and standard deviation of one or more large samples from the Binomial distribution for different values values of the parameters size and prob. Compare these to the theoretical values.

We use a large sample size (the notation 1e6 indicates 10⁶ or one million) so that the sample statistic will be very close to its theoretical value based on the parameter values. Because I work on a netbook computer with a comparatively slow processor one million is a large sample size for me. You may be able to use larger sizes.

We just saw that for an exponential distribution with parameter rate the mean is 1/rate and the standard deviation is also 1/rate. The shape of the distribution is, of course, the negative exponential shape.

The Central Limit Theorem states that for the distribution of the sample mean of i.i.d. samples of size n from a distribution with mean mu and standard deviation sigma will have mean mu and standard deviation sigma/sqrt(n) and, furthermore, the shape of the distribution will tend to a normal or Gaussian shape as n gets large.

But the result on the shape is an asymptotic result (i.e. as n gets large). Can we count on it holding for moderate values of n?

Our steps in such a simulation are:

1. Choose the values on n that we wish to examine. Because of the factor of sqrt(n) in the denominator of the standard deviation of the sample mean, we may want to use numbers that are squares, such as 1, 4, 9, 16, 25, 36, 49.
2. Choose a value of N, the number of replications. I would recommend at least 10⁵ and perhaps as much as 10⁶. The trade-offs are that a larger value allows you to get a smoother approximation to the density and more precise approximations to the theoretical quantities but large values of N take longer. Start with a moderately large value and, if your computer can handle that quickly then scale up until it seems to be taking too long.
3. Choose the distribution and parameter values and the statistic of interest. For our simulation we will use the exponential (rexp) and the default value of the parameter rate. We use the function mean to evaluate the sample mean.

If we want just one value of the sample mean from an i.i.d. sample of size n = 9 we could use

mean(rexp(9, rate = 1))

[1] 1.08991

but now we want to repeat that operation N times and save the resulting values for later analysis.

Because R is a programming language we can use some of the control structures to accomplish this. Neophyte R programmers often do this in a for loop but you have to be careful exactly how you set up the loop if you want to do this effectively.

Fortunately, there is an alternative to the for loop. The replicate function in R is a convenient way of repeating a calculation, usually involving a simulation, and collecting the results.

str(mns9 <- replicate(10^5, mean(rexp(9, rate = 1))))

 num [1:100000] 0.657 1.358 0.938 0.994 1.659 ...

So now we have generated a sample of size 100,000 of means of samples of size 9 from the standard exponential distribution.

We expect

mean(mns9)

[1] 1.000712

to be close to 1 and

sd(mns9)

[1] 0.334026

to be close to 1/3.

To examine the shape of the distribution we can use an empirical density plot

library(ggplot2)
qplot(mns9, geom="density", main="Means of samples of size 9 from an exponential dist.")

or a normal probability plot, also called a qq or quantile-quantile plot.

qplot(sample=mns9, main="Means of samples of size 9 from an exponential dist.")

That last plot takes a long time to produce because there are so many points to plot. There are techniques, to be discussed later, to thin the plot and avoid plotting so many points.