R help - Apr 2012 - automating a script to read a file

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

On Mon, Apr 23, 2012 at 04:02:45PM -0400, Steve_Friedman at nps.gov
wrote:
> 
> 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.
Hi.

If you read the parameters from a file and put them to a matrix,
then all the plots may be produced using a loop like the following.

  #some parameters
  p <- matrix(1:18, nrow=3, ncol=6)
  for (i in 1:nrow(p)) {
      plot.spdf(p[i, 1], p[i, 2], p[i, 3], p[i, 4], p[i, 5], p[i, 6])
      readline("press Enter to continue")
  }

If you use pdf() for sending the graphics to a file, then remove
the "readline" command.

Hope this helps.

Petr Savicky.

On Tue, Apr 24, 2012 at 6:00 PM, Steve Friedman <skfglades at gmail.com>
wrote:
> The original post does contain the function plot.spdf and some sample data.
That's great. And who's responsibility is it to make sure we have
context for THIS post ? This is a mailing list, not a forum and MY
mail client doesn't keep track of YOUR postings.

BTW, you are the one seeking our help, so even if you think it was a
picky request, a proper response would be to simply post your
function+data set. Not send us looking for some "original post".
> I thought I also made it clear that I understood that hist and histogram
are
> not the same.
No you didn't. But now IT IS clear you don't realize grid and base
graphics don't mix.
> I am looking for a way to use the alternatives in a function similar to but
> not identical to doubleYScale.
LatticeExtra functions expect any lattice plots, including local (user
defined) panel/plot functions, if only you would have posted plot.spdf
maybe you would have gotten a solution by now...
If plot.spdf is actually a base graphics plotter and you just stumbled
upon DoubleYscale somehow, ?par ?axis hist(...,add=T) are all
available to you.

Best
>
> Steve
>
> On Apr 24, 2012 7:45 PM, "ilai" <keren at math.montana.edu>
wrote:
>>
>> On Tue, Apr 24, 2012 at 2:09 PM, ?<Steve_Friedman at nps.gov>
wrote:
>>
>> > library(latticeExtra)
>> >> doubleYScale(hist("mydata", breaks=20, prob=T,
?xlim=c(-100, 2000),
>> > plot.spdf(x), use.style=FALSE)
>> >
>> > This does not work as doubleYScale expects histogram and
densityplot,
>> > ?and
>> > I'd like to use the plot.spdf routine in its place.
>>
>> Maybe it doesn't work because ?hist is base graphics, not the same
as
>> ?histogram in lattice. You did not provide plot.spdf or any data, so
>> this just a "place to start" type response.
>>
>> Cheers
>>
>>
>>
>> >
>> > Any Suggestions would be greatly appreciated.
>> >
>> > I'm working with
>> > R 2.15.0 (2012-03-12)
>> > Platform i386-pc-mingw32/ie86 (32-bit)
>> >
>> > Thanks
>> > Steve
>> >
>> >
>> > 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
>> >
>> >
>> >
>> > ? ? ? ? ? ? Petr Savicky
>> > ? ? ? ? ? ? <savicky at cs.cas.c
>> > ? ? ? ? ? ? z>
>> > To
>> > ? ? ? ? ? ? Sent by: ? ? ? ? ? ? ? ? ?r-help at r-project.org
>> > ? ? ? ? ? ? r-help-bounces at r-
>> > ?cc
>> > ? ? ? ? ? ? project.org
>> >
>> > Subject
>> > ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? Re: [R] automating a script
to
>> > read
>> > ? ? ? ? ? ? 04/23/2012 04:42 ? ? ? ? ?a file
>> > ? ? ? ? ? ? PM
>> >
>> >
>> >
>> >
>> >
>> >
>> >
>> >
>> >
>> > On Mon, Apr 23, 2012 at 04:02:45PM -0400, Steve_Friedman at
nps.gov wrote:
>> >>
>> >> 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.
>> >
>> > Hi.
>> >
>> > If you read the parameters from a file and put them to a matrix,
>> > then all the plots may be produced using a loop like the
following.
>> >
>> > ?#some parameters
>> > ?p <- matrix(1:18, nrow=3, ncol=6)
>> > ?for (i in 1:nrow(p)) {
>> > ? ? ?plot.spdf(p[i, 1], p[i, 2], p[i, 3], p[i, 4], p[i, 5], p[i,
6])
>> > ? ? ?readline("press Enter to continue")
>> > ?}
>> >
>> > If you use pdf() for sending the graphics to a file, then remove
>> > the "readline" command.
>> >
>> > Hope this helps.
>> >
>> > Petr Savicky.
>> >
>> > ______________________________________________
>> > R-help at 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.
>> >
>> > ______________________________________________
>> > R-help at 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.
>>
>> ______________________________________________
>> R-help at 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.