# 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)

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 # NO RESTRICTION
> sigep <- 0.01 + exp(alpha) # greater than
> 0.01
> sigeta <- 0.01 + exp(alpha) # greater than
> 0.01
> rho2 <- 0.8*sin(alpha) # between -0.8
> and 0.8
> rho1 <- lstar*abs(alpha)/(1+abs(alpha)) # between 0 and
> lstar
> delta <- 0.01 + exp(alpha) # 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:
>
>> 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
> 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