ограниченная оптимизация в R

Question

Я пытаюсь использовать http://rss.acs.unt.edu/Rdoc/library/stats/html/constrOptim.html в R для оптимизации в R с некоторыми заданными линейными ограничениями, но не могу понять, как решить проблему.

Например, мне нужно максимизировать $ f (x, y) = log (x) + \ frac {x ^ 2} {y ^ 2} $ с учетом ограничений $ g_1 (x, y) = x + y <1 $, $ g_2 (x, y) = x> 0 $ и $ g_3 (x, y) = y> 0 $. Как мне сделать это в R? Это всего лишь гипотетический пример. Не беспокойтесь о его структуре, вместо этого мне интересно узнать, как настроить это в R.

спасибо!

Answer 1

Настройка функции была тривиальной:

fr <- function(x) {      x1 <- x[1]
    x2 <- x[2]
    -(log(x1) + x1^2/x2^2)  # need negative since constrOptim is a minimization routine
}

Настройка матрицы ограничений была проблематичной из-за отсутствия большого количества документации, и я прибег к экспериментам. На странице справки написано «Возможная область определяется как ui% *% theta - ci> = 0». Итак, я проверил, и это, казалось, "работает":

> rbind(c(-1,-1),c(1,0), c(0,1) ) %*% c(0.99,0.001) -c(-1,0, 0)
      [,1]
[1,] 0.009
[2,] 0.990
[3,] 0.001

Поэтому я ставлю в ряд для каждого ограничения / границы:

constrOptim(c(0.99,0.001), fr, NULL, ui=rbind(c(-1,-1),  # the -x-y > -1
                                              c(1,0),    # the x > 0
                                              c(0,1) ),  # the y > 0
                                           ci=c(-1,0, 0)) # the thresholds

Для этой проблемы есть потенциальная трудность в том, что для всех значений x функция переходит в Inf при y -> 0. Я получаю максимум около x = .95 и y = 0, даже когда я нажимаю начальные значения в «угол», но я несколько подозреваю, что это не тот истинный максимум, который я бы предположил, был в «углу». РЕДАКТИРОВАТЬ: Следуя этому, я решил, что градиент может обеспечить дополнительное «направление», и добавил функцию градиента:

grr <- function(x) { ## Gradient of 'fr'
    x1 <- x[1]
    x2 <- x[2]
    c(-(1/x[1] + 2 * x[1]/x[2]^2),
       2 * x[1]^2 /x[2]^3 )
}

Это «направило» оптимизацию немного ближе к углу c (.999 ..., 0), вместо того, чтобы отойти от него, как это было для некоторых начальных значений. Я по-прежнему несколько разочарован тем, что процесс, похоже, «движется к утесу», когда начальные значения находятся близко к центру возможной области:

 constrOptim(c(0.99,0.001), fr, grr, ui=rbind(c(-1,-1),  # the -x-y > -1
                                               c(1,0),    # the x > 0
                                               c(0,1) ),  # the y > 0
                                            ci=c(-1,0, 0) )
$par
[1]  9.900007e-01 -3.542673e-16

$value
[1] -7.80924e+30

$counts
function gradient 
    2001       37 

$convergence
[1] 11

$message
[1] "Objective function increased at outer iteration 2"

$outer.iterations
[1] 2

$barrier.value
[1] NaN

Примечание: Ханс Вернер Борхерс опубликовал лучший пример в R-Help, в котором удалось получить значения углов, установив ограничение немного от края:

> constrOptim(c(0.25,0.25), fr, NULL, 
              ui=rbind( c(-1,-1), c(1,0),   c(0,1) ),  
              ci=c(-1, 0.0001, 0.0001)) 
$par
[1] 0.9999 0.0001