[R] Capture of iterative integration output

Sego, Landon H landon.sego at pnl.gov
Sun Aug 27 01:04:55 CEST 2006


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 at stat.math.ethz.ch 
> [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Terry W. Schulz
> Sent: Saturday, August 26, 2006 3:42 PM
> To: r-help at 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 at comcast.net
> (303) 526-1461
> 
> ______________________________________________
> R-help at 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.
>



More information about the R-help mailing list