[R] automating a script to read a file

[Steve_Friedman at nps.gov]

Hi,


The following script (which I did not develop) is used to calculate and
plot a skewed normal curve.  The script currently requires the user to
input six parameters, rather than reading these directly from a file.

I've been spinning wheels here, trying to figure out how to modify the
script to automate it.  I have four data sets, each in excess of 300
records that I need to process.

My initial thoughts were to use the  lapply and use a pdf graphic device to
capture the plots to do this, but my R programming skills are too limited
to determine how to best accomplish this.

If any one can provide assistance I would appreciate the help.


Thanks,

Steve


## Function set to find values in a skewed normal distribution

print("syntax:  plot.spdf(min, max, skewlocation, skewscale, skewshape,
skewmax, <skewtitle>)")
flush.console()

#  sample input data could be the following:
#    -100, 1000, 976.02, 230, -34, 0.7543
#         0,  500,  270, 350, -13, 0.7707
#    or any other data of similar form


erf <- function(z) {
        ## Chebyshev fitting formula for erf(z) from
        ##  Vetterling, W.T., , W.H. Press,  S.A. Teukolsky, and B.P.
Flannery. 1999.
        ##  Numerical Recipes: Example Book [C], Second Edition.
        ##  Cambridge University Press, NY. , Chapter 6-2.

      t <- 1.0/(1.0 + 0.5 * abs(z))
      ## use Horner's method
       ans <- (1 - t * exp(-z * z - 1.26551223
                + t * (1.00002368
                + t * (0.37409196
                + t * (0.09678418
                + t * (-0.18628806
                + t * (0.27886807
                + t * (-1.13520398
                + t * (1.48851587
                + t * (-0.82215223
                + t * (0.17087277)))))))))))
      if (z >= 0) return(ans) else return(-1 * ans)
}

pdf <- function(x) {
      return((1 / sqrt(2 * pi)) * exp(-1 * ((x * x) / 2)))
}

cdf <- function(x) {
      return( 0.5 * (1 + erf(x / sqrt(2))))
}

spdf <- function(x, skewlocation, skewscale, skewshape) {
      xmod <- (x - skewlocation) / skewscale
      return( 2 * pdf(xmod) * (cdf(xmod * skewshape)))
}


#### Plotting Function ################
plot.spdf <- function(xmin, xmax, skewlocation, skewscale, skewshape,
skewmax, skewtitle) {

    if(missing(skewtitle)) {
        plottitle <- "Skewed Probability Density Function"
    } else {
        plottitle <- skewtitle
      }

      skip <- (xmax - xmin) / 100.0
      xArray <- numeric(100)
      yArray <- numeric(100)

      for (i in 1:100){
            x <- xmin + i * skip
            y <- (spdf(x, skewlocation, skewscale, skewshape))/skewmax
            xArray[i] <- x
            yArray[i] <- y
      }

      plot(xArray,yArray, main=plottitle)

}





Steve Friedman Ph. D.
Ecologist  / Spatial Statistical Analyst
Everglades and Dry Tortugas National Park
950 N Krome Ave (3rd Floor)
Homestead, Florida 33034

Steve_Friedman at nps.gov
Office (305) 224 - 4282
Fax     (305) 224 - 4147