# [R] automating a script to read a file

Steve_Friedman at nps.gov Steve_Friedman at nps.gov
Mon Apr 23 22:02:45 CEST 2012

```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)