# Script runs 1,000 simulations drawing a sample of size 100 from a # Poisson distribution, and then plots histograms of the estimated # means and estimated variances. # # Figure caption: Histograms displaying distributions of X-bar and S # squared based on 1,000 randomly generated samples of size n = 100 # froma Poisson distribution with mean parameter mu = 10. In these # repeated samples, both X-bat and S squared have distributions that # are approximately normal. Both distributions are centered at 10 # (both estimators are unbiased) but the values of S squared # fluctuate much more than do the values of X-bar. nSimulations <- 1000 nSamples <- 100 lambda <- 10 # Generate all samples at once. simulationSamples <- rpois( nSimulations * nSamples, lambda) dim(simulationSamples) <- c( nSamples, nSimulations ) # Now that we have a matrix of poisson samples, with 1000 columns and # 100 rows, each column can be considered a simulation. sampleMeans <- apply(simulationSamples, 2, mean) sampleVariances <- apply(simulationSamples, 2, var) xValues<-seq(4, 16, by = 0.01) par(mfrow = c(1, 2)) par(xaxs = "i") par(yaxs = "i") par(mar = c(5, 3, 1, 3)) # Plot the sample means histogram. hist(sampleMeans, freq = FALSE, xlim = c(4, 16), ylim = c(0, 1.5), yaxt = "n", main = "", xlab = expression(bar(X)), ylab = "", cex.lab = 1.8, cex.axis = 1.6) # Add a line for the theoretical pdf of the sample mean. # This is normal with mean lambda and variance sampleMeansVariance = lambda / nSamples lines(xValues, dnorm(xValues, lambda, sqrt(sampleMeansVariance)), lwd = 2) # Plot the sample variances histogram. hist(sampleVariances, freq = FALSE, xlim = c(4, 16), ylim = c(0, 0.3), yaxt = "n", main = "", xlab = expression(S^{2}), ylab = "", cex.lab = 1.8, cex.axis = 1.6) # Compute the variance of the sample variances and plot the # theoretical density. sampleVariancesVariance <- (lambda/nSamples) + 2 * (lambda^2 / (nSamples - 1)) theoreticalPdf <- dnorm(xValues, lambda, sqrt(sampleVariancesVariance)) lines(xValues, theoreticalPdf, lwd = 2) # Print the figure as a postscript. dev.print(device = postscript, "8.2.eps", horizontal = TRUE)