Ответ 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. Я получаю max вокруг 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