CODEWITHSUNDEEP
roptimization
Student1981
Student1981

Reputation: 329

Is there anything wrong with nlminb in R?

I am trying to solve a minimization problem in R with nlminb as part of a statistical problem. However, there is something wrong when comparing the solution provided by nlminb with the plot of the function I am trying to minimize. This is the R-code of the objective function

library(cubature)
Objective_Function <- function(p0){
    F2 <- function(x){
        u.s2 <- x[1]
        u.c0 <- x[2]
        u.k0 <- x[3]
        s2 <- u.s2^(-1) - 1
        c0 <- u.c0^(-1) - 1
        k0 <- u.k0/p0
        L <- 1/2 * c0 * s2 - 1/c0 * log(1 - k0 * p0)
        A <- 1 - pnorm(L, mean = 1, sd = 1)
        A <- A * dgamma(k0, shape = 1, rate = 1)
        A <- A * dgamma(c0, shape = 1, rate = 1)
        A <- A * dgamma(s2, shape = 1, rate = 1)
        A * u.s2^(-2) * u.c0^(-2) * 1/p0
        }

    Pr <- cubature::adaptIntegrate(f = F2, 
        lowerLimit = rep(0, 3),
        upperLimit = rep(1, 3))$integral

    A <- 30 * Pr * (p0 - 0.1)
    B <- 30 * Pr * (1 - Pr) * (p0 - 0.1)^2
    0.4 * B + (1 - 0.4) * (-A)
    }

Following the R-command

curve(Objective_Function, 0.1, 4)

one observes a critical point close to 2. However, when one executes

nlminb(start = runif(1, min = 0.1, max = 4), 
    objective = Objective_Function,
    lower = 0.1, upper = 4)$par

the minimum of the function takes place at the point 0.6755844.

I was wondering if you could tell me where my mistake is, please. Is there any reliable R-command to solve optimization problems?

If this is a very basic question, I apologize.

Thank you for your help.

Upvotes: 0

Views: 607

Answers (1)

Dimitris Rizopoulos
Dimitris Rizopoulos

Reputation: 350

The problem is not nlminb() but the fact that you have not provided a vectorized function in curve(). You can get the correct figure using the following code, from which you see that nlminb() indeed finds the minimum:

min_par <- nlminb(start = runif(1, min = 0.1, max = 4), 
                  objective = Objective_Function,
                  lower = 0.1, upper = 4)$par

vec_Objective_Function <- function (x) sapply(x, Objective_Function)
curve(vec_Objective_Function, 0.1, 4)
abline(v = min_par, lty = 2, col = 2)

In addition, for univariate optimization you can also use function optimize(), i.e.,

optimize(Objective_Function, c(0.1, 4))

Upvotes: 1

Related Questions