Re: [R] Problem with numerical integration and optimization with BFGS

From: Prof Brian Ripley <ripley_at_stats.ox.ac.uk>
Date: Fri, 25 May 2007 04:22:27 +0100 (BST)

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

-- 
Brian D. Ripley,                  ripley_at_stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595

______________________________________________
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.
Received on Fri 25 May 2007 - 03:31:20 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Fri 25 May 2007 - 06:31:35 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.