R help - Aug 2006 - 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.3989422804014327/
                           (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

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.
>