From: Sego, Landon H <landon.sego_at_pnl.gov>

Date: Sun 27 Aug 2006 - 09:04:55 EST

R-help@stat.math.ethz.ch 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. Received on Sun Aug 27 09:08:23 2006

Date: Sun 27 Aug 2006 - 09:04:55 EST

capture <- double(n) # creates a vector of length n, all with values of 0

capture[i] <- integrate( _your arguments_ )$value

"capture" should be the same as your "z".

> -----Original Message-----

*> From: r-help-bounces@stat.math.ethz.ch
**> [mailto:r-help-bounces@stat.math.ethz.ch] On Behalf Of Terry W. Schulz
**> Sent: Saturday, August 26, 2006 3:42 PM
**> To: r-help@stat.math.ethz.ch
**> Subject: [R] Capture of iterative integration output
**>
**> Hello,
**>
**> I am a novice R user and am having difficulty retrieving the
**> values from 21 iterations of the R function integrate.
**> The only way I have found is to do a write.table and then a
**> read.table as shown in the code below. I would rather
**> capture the 21 values inside the braces ( sapply might work,
**> but I can't set it up without getting an error in function)
**> so I could compute the negative log-likelihood inside the
**> braces. My goal is to then use optim to find the best
**> fitting mu and sigma. In the code that follows mu and sigma
**> are given, but normally they would not be known.
**>
**> Any help on y$value capture would be appreciated. I would be
**> ecstatic for help on implementation of the optim function
**> using finite differences for the following code.
**>
**> By the way in case anyone is interested function(x) is the
**> Poisson-lognormal pdf.
**>
**> q <- read.table(file="c:/Bayes/FigA3-1.prn",header=TRUE)
**> attach(q)
**> V <- c(volume)
**> c <- c(count)
**> n <- length(count)
**> mu <- -8.202900475
**> sigma <- 1.609437912
**> for(i in c(1:n))
**> {
**> integrand <- function(x)
**> {(x*V[[i]])^c[[i]]*exp(-x*V[[i]])/factorial(c[[i]])*0.39894228
*

04014327/

> (x*sigma)*exp(-.5*((log(x)-mu)/sigma)^2)}

*> y <- integrate(integrand, lower = 0, upper = Inf)
**> write.table(y[[1]],file="c:/Bayes/PLNA3-1.out",append=TRUE,
**> quote=FALSE,row.names=FALSE,col.names=FALSE)
**> }
**> z <-read.table(file="c:/Bayes/PLNA3-1.out")
**> z <- c(z)
**> LL <- log(c(z$V1))
**> #negative Log-Likelihood
**> print(-sum(LL),digits=6)
**> file.remove("c:/Bayes/PLNA3-1.out")
**>
**> --
**> Terry W. Schulz
**> Applied Statistician
**> 1218 Pomegranate Lane
**> Golden, CO 80401
**> USA
**>
**> terrywschulz@comcast.net
**> (303) 526-1461
**>
**> ______________________________________________
**> R-help@stat.math.ethz.ch 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@stat.math.ethz.ch 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. Received on Sun Aug 27 09:08:23 2006

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.1.8, at Sun 27 Aug 2006 - 12:23:43 EST.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help.
Please read the posting
guide before posting to the list.
*