# R macros for Math 324, originally written spring 2005.

rdist <- function(num.stat, stat, dist, ...) {
  # num.stat is the number of statistics `stat' to be generated.
  # Some choices for `stat':  mean, median, sum, sd, var, min, max
  # Some choices for `dist':  
  #    rnorm, rexp, rcauchy, rlaplace, rt, rchisq, rf
  # ...: optional arguments to `dist'
  # The first argument (`n') of `dist' is the number of observations
  #    used to compute each stat, and is necessary.
  x <- NA
  for (i in 1:num.stat) { x <- c(x, match.fun(stat)(match.fun(dist)(...))) } 
  x[2:length(x)]
}


plot.dist <- function(xmin, xmax, dist, ...) {
  # `xmin' and `xmax' represent the limiting values on the x-axis.
  # When plotting the pdf, some choices for `dist' are
  #    dnorm, dexp, dcauchy, dlaplace, dt, dchisq, df, dunif, dbinom, dgeo, dpois
  # When plotting the cdf, some choices for `dist' are
  #    pnorm, pexp, pcauchy, plaplace, pt, pchisq, pf, punif, pbinom, pgeo, ppois
  # ...: optional arguments to `dist', EXCLUDING the first argument.
  num.grid.points <- 10000;    options(warn= -1)
  x <- xmin + (xmax-xmin)*(0:num.grid.points)/num.grid.points
  ymax <- max( match.fun(dist)(x, ...) )
  discrete.pdf <- sum( match.fun(dist)(x[floor(x)!=ceiling(x)], ...) )==0
  if (discrete.pdf) { x <- ceiling(xmin):floor(xmax)
     plot(x, (match.fun(dist)(x, ...)), ylim=c(0,ymax), type="h", ylab="")  }
  else { plot(x, (match.fun(dist)(x, ...)), ylim=c(0,ymax), type="s", ylab="") }
}
                                                                                                   
trunc.hist <- function(xmin, xmax, x, ...) {
  # Produces a truncated histogram.
  # When constructing the histogram, exclude values
  #    outside the domain from `xmin' to `xmax'.
  # `x' is the original data set.
  # ...: optional arguments to `hist'.
  y <- x[ x>xmin & x<xmax ]
  hist( y, ...)
}
  


rlaplace <- function(n, sd=1) {
  # Random generation for the Laplace (double exponential) distribution 
  #   with standard devation `sd'.
  # `n' is the number of observations.
  rexp(n, (1/sd)) * sample(c(-1, 1), n, replace=T) / sqrt(2)
} 

dlaplace <- function(x, sd=1) {
  # Probability density function (pdf) for the Laplace 
  #   (double exponential) distribution with standard devation `sd'.
  exp(-abs(x)*sqrt(2)/sd)/(sd*sqrt(2))
}

plaplace <- function(x, sd=1) {
  # Cumulative distribution function (cdf) for the Laplace 
  #   (double exponential) distribution with standard devation `sd'.
  (x>0) - ((x>0)*2-1) * exp(-abs(x)*sqrt(2)/sd)/2 
}