Re: [R] Capture of iterative integration output

From: Sego, Landon H <landon.sego_at_pnl.gov>
Date: Sun 27 Aug 2006 - 09:04:55 EST


Here's one way to do this:

Before you begin the loop, define a vector to capture the output. Then each iteration of the loop will be assigned to the corresponding element in that vector:

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

for (i in 1:n) {

  _ some code _

  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.