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_fried...@nps.gov Office (305) 224 - 4282 Fax (305) 224 - 4147 ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.