# S-PLUS file for Writing Assignment #1: Understanding Random Sampling
#
# Save this file as a .ssc file (for S-PLUS script),
# either on the desktop (in Windows) or to a floppy disk.
#
# You may load these functions into S-PLUS by doing the following:
#
# 1. Start S-PLUS
# 2. If you saved this on a floppy, in the command window,
#    after the > prompt, type
#
#    > source("a:filename.ssc")
#
#    where filename is the name of the file.
#
#    If it is on the desktop, choose Open from the File menu, look in the desktop,
#    and open the file.  Then, from the Script menu, choose run.
#
# After running the script, you may use these functions to take many samples
# of any size desired from the five benchmark populations, find the sample means
# of the samples, and summarize the data.
#
# The names of the functions for sampling from the benchmark populations are:
# 
# 1. snorm
# 2. sskew
# 3. sbimodal
# 4. sdsym
# 5. sdskew
#
# For example, to take a sample of size 1000 from the skewed population
# and save the sample in an object called x, type
#
# > x <- sskew(1000)
#
# You could find the mean of this sample by typing
#
# > mean(x)
#
# You could plot a Q-Q plot versus the normal distribution by typing
#
# > qqnorm(x)
#
# If you want to take 10,000 separate samples of size 25 from the bimodal population,
# find the mean of each sample, and save these 10,000 sample means
# in an object called x, type
#
# > x <- get.means(bimodal,10000,25)
#
# If you then wanted to know how many of these sample means were between
# 98 and 102, you could type
#
# > sum( (x > 98) & (x < 102) )

snorm <- function(n)
{ 
  return(rnorm(n,100,10))
}

sskew <- function(n)
{
  return(rexp(n,rate=0.1) + 90)
}

sbimodal <- function(n)
{
  b <- rbinom(n,1,0.5)
  x1 <- rnorm(n,91,sqrt(19))
  x2 <- rnorm(n,109,sqrt(19))
  return(b*x1 + (1-b)*x2)
}

sdsym <- function(n)
{
 return(rbinom(n,1,.5) * 20 + 90)
}

sdskew <- function(n)
{
  return(rbinom(n,1,1/101) * 101 + 99)
}

get.means <- function(f,nsamples,size=25)
{
  ans <- matrix(NA,nsamples,size)
  for( i in 1:nsamples )
  {
    ans[i,] <- f(size)
  }
  return( apply(ans,1,mean) )
}