[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)
Homestead, Florida 33034
Steve_Friedman at nps.gov
Office (305) 224 - 4282
Fax (305) 224 - 4147
More information about the R-help
mailing list