[R] minimizing a numerical integration

Hans W Borchers hwborchers at googlemail.com
Sat Feb 16 09:12:44 CET 2013


> Dear all,
> 
> I am a new user to R and I am using pracma and nloptr libraries to minimize
> a numerical integration subject to a single constraint . The integrand
> itself is somehow a complicated function of x and y that is computed
> through several steps.  i formulated the integrand in a separate function
> called f which is a function of x &y. I want to find the optimal value of x
> such that the integration over y is minimum. Here is my code, it is not
> working. I got an error: 'y' is missing
> 
> library('nloptr')
>  library('pracma')

Guess there's no need to use pracma or nloptr. `eval_f0` returns the error
message as function `f(x, y)` when called in an integration routine cannot
decide which variable to integrate over.

You don't provide code, so here is a simple example:

    f <- function(x, y) sin(y) * cos(x * y)
    eval_f0 <- function(x) integrate(function(y) f(x, y), 0, pi)$value

    optimize(eval_f0, c(0, 2*pi))
    ## minimum: 1.652188
    ## objective: -0.844125

In your code, x is a scalar. But if x is a vector, applying nloptr might be a
choice:

    f <- function(x, y) sin(y) * cos((x[1]+x[2])*y)

    eval_f0 <- function(x) integrate(function(y) f(x, y), 0, pi)$value
    eval_g0 <- function(x) x[1]^2 + x[2]^2 - 1  # i.e., sum(x^2) <= 1

    nloptr( x0=c(0.5, 0.5),
            eval_f=eval_f0,
            lb = c(0, 0),
            ub = c(1, 1),
            eval_g_ineq = eval_g0,
            opts = list("algorithm"="NLOPT_LN_COBYLA", "maxeval"=1000))
    ## Optimal value of objective function:  -0.733744658465974 
    ## Optimal value of controls: 0.707091 0.7071225
    
Hans Werner

> f <- function(x,y) {#here i should put the commands representing my function
> return(   )
> }
> 
>  #constraint function
> eval_g0 <- function(x) {
>     return( )
> }
> 
> # objective function
> eval_f0 <- function(x) {
> romberg(f, 0.5, 0.5001)}
> 
> ARL1 <- nloptr( x0=c(0.653),
>                 eval_f=eval_f0,
>                 lb = c(0),
>                 ub = c(6),
>                 eval_g_ineq = eval_g0,
>                 opts = list("algorithm"="NLOPT_LN_COBYLA", "maxeval"=1000),
> 
>  Thanks



More information about the R-help mailing list