[IPSUR-commits] r118 - in pkg/IPSUR: R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Jan 4 15:19:50 CET 2010
Author: gkerns
Date: 2010-01-04 15:19:50 +0100 (Mon, 04 Jan 2010)
New Revision: 118
Added:
pkg/IPSUR/R/clt.R
pkg/IPSUR/man/clt.Rd
Log:
added central limit theorem exercises
Added: pkg/IPSUR/R/clt.R
===================================================================
--- pkg/IPSUR/R/clt.R (rev 0)
+++ pkg/IPSUR/R/clt.R 2010-01-04 14:19:50 UTC (rev 118)
@@ -0,0 +1,311 @@
+##################################################
+# The Central Limit Theorem
+# want to investigate how the distribution of
+# x-bar changes as the sample size gets large
+
+population <- "rt" # pop'n distribution is Student's t
+r <- 3 # degrees of freedom parameter
+
+sample.size <- 2 # sample size
+
+N.iter <- 100000 # number of simulated xbar's
+
+
+clt1 <- function(population = "rt",
+ r = 3,
+ sample.size = 2,
+ N.iter = 100000){
+
+#################################################################
+# initialize variables
+population <- get(population, mode = "function")
+xbar <- rep(0, N.iter)
+graphics.off()
+
+curve( dt(x, df = r ),
+ xlim = c(-5,5),
+ xlab = "Support Set",
+ ylab = "Density",
+ lwd = 2,
+ main = "The Population Distribution \n (while we're waiting)" )
+abline( h = 0 , col = "grey" )
+
+
+########################################
+# Label the plot with mu
+text( 5,
+ dt(0, df = r )*0.9,
+ bquote( mu ==.(0) ),
+ cex = 1.5,
+ pos = 2 )
+
+# Label the plot with sigma^2
+text( 5,
+ dt(0, df = r )*0.8,
+ bquote( sigma^2 ==.(r/(r-2)) ),
+ cex = 1.5,
+ pos = 2 )
+
+
+#################################################
+# simulate xbar's
+
+xbar <- rowMeans( matrix(population(sample.size * N.iter, df = r),
+ nrow = N.iter)
+ )
+
+# Find mean and variance of xbar
+xbar.mean <- round( mean( xbar ), 4)
+xbar.var <- round( var( xbar ), 4)
+
+# window for graph
+low <- xbar.mean - 3*sqrt(xbar.var)
+up <- xbar.mean + 3*sqrt(xbar.var)
+
+dev.new()
+dev.set(3)
+# Draw histogram of simulated x-bars
+hist( xbar,
+ breaks = 280,
+ xlim = c(low,up),
+ xlab = "",
+ prob = TRUE,
+ main = "Sampling Distribution of X-bar",
+ sub = "Click to see Limiting Normal Density (in red)")
+
+########################################
+# Label the histogram with mean(xbar)
+text( up,
+ dnorm(xbar.mean, mean = xbar.mean, sd = sd(xbar)),
+ bquote( mean(xbar)==.(xbar.mean) ),
+ cex = 1,
+ pos = 2 )
+
+# Label the histogram with var(xbar)
+text( up,
+ dnorm(xbar.mean, mean = xbar.mean, sd = sd(xbar))*0.9,
+ bquote( var(xbar) ==.(xbar.var) ),
+ cex = 1,
+ pos = 2 )
+
+# Label the histogram with n
+text( up,
+ dnorm(xbar.mean, mean = xbar.mean, sd = sd(xbar))*0.8,
+ bquote( n ==.(sample.size) ),
+ cex = 1.85,
+ pos = 2 )
+
+######################################
+# Draw limiting Normal curve
+z <- locator( n = 1 )
+curve( dnorm(x, mean = xbar.mean, sd = sd(xbar)),
+ lwd = 2,
+ col = "red",
+ add = T )
+}
+
+
+
+
+
+
+
+clt2 <- function(population = "runif",
+ a = 0,
+ b = 10,
+ sample.size = 2,
+ N.iter = 100000){
+
+#################################################################
+# initialize variables
+population <- get(population, mode = "function")
+xbar <- rep(0, N.iter)
+graphics.off()
+
+curve( dunif(x, min = a, max = b ),
+ xlim = c(a-1,b+1), ylim = c(0, 1.3/(b-a)),
+ xlab = "Support Set",
+ ylab = "Density",
+ lwd = 2,
+ main = "The Population Distribution \n (while we're waiting)" )
+abline( h = 0 , col = "grey" )
+
+
+########################################
+# Label the plot with mu
+text( (a+b)/2,
+ 0.9/(b-a),
+ bquote( mu ==.((a+b)/2) ),
+ cex = 1.5,
+ pos = 1 )
+
+# Label the plot with sigma^2
+text( (a+b)/2,
+ 0.8/(b-a),
+ bquote( sigma^2 ==.( (b-a)^2/12 ) ),
+ cex = 1.5,
+ pos = 1 )
+
+
+#############################################
+# simulate xbar's
+xbar <-rowMeans(matrix(population(sample.size * N.iter, min = a, max = b),
+ nrow = N.iter)
+ )
+
+# Find mean and variance of xbar
+xbar.mean <- round( mean( xbar ), 4)
+xbar.var <- round( var( xbar ), 4)
+
+# window for graph
+low <- xbar.mean - 3*sqrt(xbar.var)
+up <- xbar.mean + 3*sqrt(xbar.var)
+
+dev.new()
+dev.set(3)
+# Draw histogram of simulated x-bars
+hist( xbar,
+ breaks = 80,
+ xlim = c(low,up),
+ xlab = "",
+ prob = TRUE,
+ main = "Sampling Distribution of X-bar",
+ sub = "Click to see Limiting Normal Density (in red)")
+
+########################################
+# Label the histogram with mean(xbar)
+text( up,
+ dnorm(xbar.mean, mean = xbar.mean, sd = sd(xbar)),
+ bquote( mean(xbar)==.(xbar.mean) ),
+ cex = 1,
+ pos = 2 )
+
+# Label the histogram with var(xbar)
+text( up,
+ dnorm(xbar.mean, mean = xbar.mean, sd = sd(xbar))*0.9,
+ bquote( var(xbar)==.(xbar.var) ),
+ cex = 1,
+ pos = 2 )
+
+# Label the histogram with n
+text( up,
+ dnorm(xbar.mean, mean = xbar.mean, sd = sd(xbar))*0.8,
+ bquote( n ==.(sample.size) ),
+ cex = 1.85,
+ pos = 2 )
+
+######################################
+# Draw limiting Normal curve
+z <- locator( n = 1 )
+curve( dnorm(x, mean = xbar.mean, sd = sd(xbar)),
+ lwd = 2,
+ col = "red",
+ add = T )
+}
+
+
+clt3 <- function(population = "rgamma",
+ alpha = 1.21,
+ theta = 2.37,
+ sample.size = 2,
+ N.iter = 100000){
+
+#################################################################
+# initialize variables
+population <- get(population, mode = "function")
+xbar <- rep(0, N.iter)
+graphics.off()
+
+
+curve( dgamma(x, shape = alpha, scale = theta ),
+ xlim = c(0, alpha*theta*(1 + 3*theta)),
+ xlab = "Support Set",
+ ylab = "Density",
+ lwd = 2,
+ main = "The Population Distribution \n (while we're waiting)" )
+abline( h = 0 , col = "grey" )
+
+f = function(x){dgamma(x, shape = alpha, scale = theta )}
+
+OPT = optimize( f,
+ interval = c(0, alpha*theta*(1 + 3*theta)),
+ maximum = TRUE)
+
+########################################
+# Label the plot with mu
+text( alpha*theta*(1 + 2*theta),
+ (OPT$objective)*0.9,
+ bquote( mu ==.(alpha*theta )),
+ cex = 1.5,
+ pos = 1 )
+
+# Label the plot with sigma^2
+text( alpha*theta*(1 + 2*theta),
+ (OPT$objective)*0.8,
+ bquote( sigma^2 ==.( alpha*theta^2 ) ),
+ cex = 1.5,
+ pos = 1 )
+
+
+#############################################
+# simulate xbar's
+xbar <- rowMeans(matrix(population(sample.size * N.iter, shape = alpha, scale = theta),
+ nrow = N.iter)
+ )
+
+# Find mean and variance of xbar
+xbar.mean <- round( mean( xbar ), 4)
+xbar.var <- round( var( xbar ), 4)
+
+# window for graph
+low <- xbar.mean - 3*sqrt(xbar.var)
+up <- xbar.mean + 3*sqrt(xbar.var)
+
+dev.new()
+dev.set(3)
+# Draw histogram of simulated x-bars
+hist( xbar,
+ breaks = 80,
+ xlim = c(low,up),
+ xlab = "",
+ prob = TRUE,
+ main = "Sampling Distribution of X-bar",
+ sub = "Click to see Limiting Normal Density (in red)")
+
+########################################
+# Label the histogram with mean(xbar)
+text( up,
+ dnorm(xbar.mean, mean = xbar.mean, sd = sd(xbar)),
+ bquote( mean(xbar)==.(xbar.mean) ),
+ cex = 1,
+ pos = 2 )
+
+# Label the histogram with var(xbar)
+text( up,
+ dnorm(xbar.mean, mean = xbar.mean, sd = sd(xbar))*0.9,
+ bquote( var(xbar)==.(xbar.var) ),
+ cex = 1,
+ pos = 2 )
+
+# Label the histogram with n
+text( up,
+ dnorm(xbar.mean, mean = xbar.mean, sd = sd(xbar))*0.8,
+ bquote( n ==.(sample.size) ),
+ cex = 1.85,
+ pos = 2 )
+
+######################################
+# Draw limiting Normal curve
+z = locator( n = 1 )
+curve( dnorm(x, mean = xbar.mean, sd = sd(xbar)),
+ lwd = 2,
+ col = "red",
+ add = T )
+
+}
+
+
+
+
+
+
Added: pkg/IPSUR/man/clt.Rd
===================================================================
--- pkg/IPSUR/man/clt.Rd (rev 0)
+++ pkg/IPSUR/man/clt.Rd 2010-01-04 14:19:50 UTC (rev 118)
@@ -0,0 +1,38 @@
+\name{The Central Limit Theorem}
+\alias{The Central Limit Theorem}
+\alias{clt1}
+\alias{clt2}
+\alias{clt3}
+
+\title{Investigating the Central Limit Theorem}
+\description{
+ These functions were written for students to investigate the Central Limit Theorem. For more information, see the exercises at the end of the chapter "Sampling Distributions" in IPSUR.
+}
+
+\usage{
+clt1(population = "rt", r = 3, sample.size = 2, N.iter = 100000)
+clt2(population = "runif", a = 0, b = 10, sample.size = 2, N.iter = 100000)
+clt3(population = "rgamma", alpha = 1.21, theta = 2.37, sample.size = 2, N.iter = 100000)
+}
+
+\arguments{
+ \item{population}{the name of a population distribution, in its random generator form.}
+ \item{sample.size}{the sample size.}
+ \item{N.iter}{the number of samples desired.}
+ \item{r}{the degrees of freedom for Student's t distribution.}
+ \item{a}{the minimum value of a continuous uniform distribution.}
+ \item{b}{the maximum value of a continuous uniform distribution.}
+ \item{alpha}{the shape parameter of a gamma distribution.}
+ \item{theta}{the scale parameter of a gamma distribution.}
+}
+
+\details{
+ When the functions are called a plot window opens to show a graph of the PDF of the population distribution. On the display are shown numerical values of the population mean and variance. The computer simulates random samples of size \code{sample.size} from the distribution a total of
+\code{N.iter} times, and sample means are calculated for each sample. Next follows a histogram of the simulated sample means, which closely approximates the sampling distribution of the sample mean.
+Also shown are the sample mean and sample variance of all of the simulated sample means. As a final step, when the user clicks the second plot, a normal curve with the same mean and variance as the simulated sample means is superimposed over the histogram.
+}
+
+
+\author{G. Jay Kerns \email{gkerns at ysu.edu}}
+
+\keyword{misc}
More information about the IPSUR-commits
mailing list