[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