[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