[R] Problem with numerical integration and optimization with BFGS
Ravi Varadhan
rvaradhan at jhmi.edu
Fri May 25 15:56:48 CEST 2007
Deepankar,
If the problem seems to be in the evaluation of numerical quadrature part,
you might want to try quadrature methods that are better suited to
integrands with strong peaks. The traditional Gaussian quadrature methods,
even their adaptive versions such as Gauss-Kronrod, are not best suited for
integrating because they do not explicitly account for the "peakedness" of
the integrand, and hence can be inefficient and inaccurate. See the article
below:
http://citeseer.ist.psu.edu/cache/papers/cs/18996/http:zSzzSzwww.sci.wsu.edu
zSzmathzSzfacultyzSzgenzzSzpaperszSzmvn.pdf/genz92numerical.pdf
Alan Genz has worked on this problem a lot and has a number of computational
tools available. I used some of them when I was working on computing Bayes
factors for binomial regression models with different link functions. If
you are interested, check the following:
http://www.math.wsu.edu/faculty/genz/software/software.html.
For your immediate needs, there is an R package called "mnormt" that has a
function for computing integrals under a multivariate normal (and
multivariate t) densities, which is actually based on Genz's Fortran
routines. You could try that.
Ravi.
----------------------------------------------------------------------------
-------
Ravi Varadhan, Ph.D.
Assistant Professor, The Center on Aging and Health
Division of Geriatric Medicine and Gerontology
Johns Hopkins University
Ph: (410) 502-2619
Fax: (410) 614-9625
Email: rvaradhan at jhmi.edu
Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
----------------------------------------------------------------------------
--------
-----Original Message-----
From: r-help-bounces at stat.math.ethz.ch
[mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Deepankar Basu
Sent: Friday, May 25, 2007 12:02 AM
To: Prof Brian Ripley
Cc: r-help at stat.math.ethz.ch
Subject: Re: [R] Problem with numerical integration and optimization with
BFGS
Prof. Ripley,
The code that I provided with my question of course does not contain
code for the derivatives; but I am supplying analytical derivatives in
my full program. I did not include that code with my question because
that would have added about 200 more lines of code without adding any
new information relevant for my question. The problem that I had pointed
to occurs whether I provide analytical derivatives or not to the
optimization routine. And the problem was that when I use the "BFGS"
method in optim, I get an error message saying that the integrals are
probably divergent; I know, on the other hand, that the integrals are
convergent. The same problem does not arise when I instead use the
Nelder-Mead method in optim.
Your suggestion that the expression can be analytically integrated
(which will involve "pnorm") might be correct though I do not see how to
do that. The integrands are the bivariate normal density functions with
one variable replaced by known quantities while I integrate over the
second.
For instance, the first integral is as follows: the integrand is the
bivariate normal density function (with general covariance matrix) where
the second variable has been replaced by
y[i] - rho1*y[i-1] + delta
and I integrate over the first variable; the range of integration is
lower=-y[i]+rho1*y[i-1]
upper=y[i]-rho1*y[i-1]
The other two integrals are very similar. It would be of great help if
you could point out how to integrate the expressions analytically using
"pnorm".
Thanks.
Deepankar
On Fri, 2007-05-25 at 04:22 +0100, Prof Brian Ripley wrote:
> You are trying to use a derivative-based optimization method without
> supplying derivatives. This will use numerical approoximations to the
> derivatives, and your objective function will not be suitable as it is
> internally using adaptive numerical quadrature and hence is probably not
> close enough to a differentiable function (it may well have steps).
>
> I believe you can integrate analytically (the answer will involve pnorm),
> and that you can also find analytical derivatives.
>
> Using (each of) numerical optimization and integration is a craft, and it
> seems you need to know more about it. The references on ?optim are too
> advanced I guess, so you could start with Chapter 16 of MASS and its
> references.
>
> On Thu, 24 May 2007, Deepankar Basu wrote:
>
> > Hi R users,
> >
> > I have a couple of questions about some problems that I am facing with
> > regard to numerical integration and optimization of likelihood
> > functions. Let me provide a little background information: I am trying
> > to do maximum likelihood estimation of an econometric model that I have
> > developed recently. I estimate the parameters of the model using the
> > monthly US unemployment rate series obtained from the Federal Reserve
> > Bank of St. Louis. (The data is freely available from their web-based
> > database called FRED-II).
> >
> > For my model, the likelihood function for each observation is the sum of
> > three integrals. The integrand in each of these integrals is of the
> > following form:
> >
> > A*exp(B+C*x-D*x^2)
> >
> > where A, B, C and D are constants, exp() is the exponential function and
> > x is the variable of integration. The constants A and D are always
> > positive; B is always negative, while there is no a priori knowledge
> > about the sign of C. All the constants are finite.
> >
> > Of the three integrals, one has finite limits while the other two have
> > the following limits:
> >
> > lower = -Inf
> > upper = some finite number (details can be found in the code below)
>
> Try integrating that analytically by change of variable to a normal CDF.
>
>
> > My problem is the following: when I try to maximize the log-likelihood
> > function using "optim" with method "BFGS", I get the following error
> > message (about the second integral):
> >
> >> out <- optim(alpha.start, LLK, gr=NULL, method="BFGS", y=urate$y)
> > Error in integrate(f3, lower = -Inf, upper = upr2) :
> > the integral is probably divergent
> >
> > Since I know that all the three integrals are convergent, I do not
> > understand why I am getting this error message. My first question: can
> > someone explain what mistake I am making?
> >
> > What is even more intriguing is that when I use the default method
> > (Nelder-Mead) in "optim" instead of BFGS, I do not get any such error
> > message. Since both methods (Nelder-Mead and BFGS) will need to evaluate
> > the integrals, my second question is: why this difference?
> >
> > Below, I am providing the code that I use. Any help will be greatly
> > appreciated.
> >
> >
> > Deepankar
> >
> >
> > ************ CODE START *******************
> >
> >
> >
> > #############################
> > # COMPUTING THE LOGLIKELIHOOD
> > # USING NUMERICAL INTEGRALS
> > #############################
> >
> > LLK <- function(alpha, y) {
> >
> > n <- length(y)
> > lglik <- numeric(n) # TO BE SUMMED LATER TO GET THE LOGLIKELIHOOD
> >
> > lambda <- numeric(n-1) # GENERATING *lstar*
> > for (i in 1:(n-1)) { # TO USE IN THE
> > lambda[i] <- y[i+1]/y[i] # RE-PARAMETRIZATION BELOW
> > }
> > lstar <- (min(lambda)-0.01)
> >
> >
> > # NOTE RE-PARAMETRIZATION
> > # THESE RESTRICTIONS EMERGE FROM THE MODEL
> >
> > muep <- alpha[1] # NO RESTRICTION
> > sigep <- 0.01 + exp(alpha[2]) # greater than
> > 0.01
> > sigeta <- 0.01 + exp(alpha[3]) # greater than
> > 0.01
> > rho2 <- 0.8*sin(alpha[4]) # between -0.8
> > and 0.8
> > rho1 <- lstar*abs(alpha[5])/(1+abs(alpha[5])) # between 0 and
> > lstar
> > delta <- 0.01 + exp(alpha[6]) # greater than
> > 0.01
> >
> >
> > ##########################################
> > # THE THREE FUNCTIONS TO INTEGRATE
> > # FOR COMPUTING THE LOGLIKELIHOOD
> > ##########################################
> >
> > denom <- 2*pi*sigep*sigeta*(sqrt(1-rho2^2)) # A CONSTANT TO BE USED
> > # FOR DEFINING THE
> > # THREE FUNCTIONS
> >
> >
> > f1 <- function(z1) { # FIRST FUNCTION
> >
> > b11 <- ((z1-muep)^2)/((-2)*(1-rho2^2)*(sigep^2))
> > b12 <-
> > (rho2*(z1-muep)*(y[i]-y[i-1]+delta))/((1-rho2^2)*sigep*sigeta)
> > b13 <- ((y[i]-y[i-1]+delta)^2)/((-2)*(1-rho2^2)*(sigeta^2))
> >
> > return((exp(b11+b12+b13))/denom)
> > }
> >
> > f3 <- function(z3) { # SECOND FUNCTION
> >
> > b31 <- ((y[i]-rho1*y[i-1]-muep)^2)/((-2)*(1-rho2^2)*(sigep^2))
> > b32 <-
> > (rho2*(y[i]-rho1*y[i-1]-muep)*z3)/((1-rho2^2)*sigep*sigeta)
> > b33 <- ((z3)^2)/((-2)*(1-rho2^2)*(sigeta^2))
> >
> > return((exp(b31+b32+b33))/denom)
> > }
> >
> > f5 <- function(z5) { # THIRD FUNCTION
> >
> > b51 <- ((-y[i]+rho1*y[i-1]-muep)^2)/((-2)*(1-rho2^2)*sigep^2)
> > b52 <-
> > (rho2*(-y[i]+rho1*y[i-1]-muep)*(z5))/((1-rho2^2)*sigep*sigeta)
> > b53 <- ((z5)^2)/((-2)*(1-rho2^2)*(sigeta^2))
> >
> > return((exp(b51+b52+b53))/denom)
> > }
> >
> >
> > for (i in 2:n) { # START FOR LOOP
> >
> > upr1 <- (y[i]-rho1*y[i-1])
> > upr2 <- (y[i]-y[i-1]+delta)
> >
> > # INTEGRATING THE THREE FUNCTIONS
> > out1 <- integrate(f1, lower = (-1)*upr1, upper = upr1)
> > out3 <- integrate(f3, lower = -Inf, upper = upr2)
> > out5 <- integrate(f5, lower= -Inf, upper = upr2)
> >
> > pdf <- (out1$val + out3$val + out5$val)
> >
> > lglik[i] <- log(pdf) # LOGLIKELIHOOD FOR OBSERVATION i
> >
> > } # END FOR LOOP
> >
> > return(-sum(lglik)) # RETURNING NEGATIVE OF THE LOGLIKELIHOOD
> > # BECAUSE optim DOES MINIMIZATION BY DEFAULT
> > }
> >
> > ***************** CODE ENDS *********************************
> >
> > Then I use:
> >
> >> urate <- read.table("~/Desktop/UNRATE1.txt", header=TRUE) # DATA
> >> alpha.start <- c(0.5, -1, -1, 0, 1, -1) # STARTING VALUES
> >> out <- optim(alpha.start, LLK, gr=NULL, y=urate$y) # THIS GIVES NO
> > ERROR
> >
> > or
> >
> >> out <- optim(alpha.start, LLK, gr=NULL, method="BFGS", y=urate$y)
> > Error in integrate(f3, lower = -Inf, upper = upr2) :
> > the integral is probably divergent
> >
> > ______________________________________________
> > 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.
> >
>
______________________________________________
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