Оптимизация с ограничениями

Question

Я работаю с выходом модели, в которой есть оценки параметров, которые не могут следовать априорным ожиданиям. Я хотел бы написать функцию, которая заставляет эти оценки полезности соответствовать этим ожиданиям. Для этого функция должна минимизировать сумму квадрата отклонения между начальными значениями и новыми оценками. Так как у нас априори ожидания, оптимизация должна быть с учетом следующих ограничений:

B0 < B1 
B1 < B2 
... 
Bj < Bj+1 

Например, оценки параметров сырья ниже flipflopped для B2 и B3. Столбцы Delta и Delta^2 показывают отклонение между исходной оценкой параметров и новым коэффициентом. Я пытаюсь свести к минимуму столбец Delta^2. Я закодирован это в Excel и показано, как Solver в Excel позволит оптимизировать эту проблему, обеспечивающую набор ограничений:

Beta BetaRaw Delta Delta^2 BetaNew 
B0  1.2  0  0   1.2 
B1  1.3  0  0   1.3 
B2  1.6  -0.2  0.04  1.4 
B3  1.4  0  0   1.4 
B4  2.2  0  0   2.2 

После прочтения ?optim и ?constrOptim, я не в состоянии, как обращал внимания установить это в Р. Я уверен, что я просто немного плотный, но могу использовать некоторые указатели в правильном направлении!

4/24/2012 - Добавлена ​​награда, так как я недостаточно умен, чтобы перевести первый ответ.

Вот код R, который должен быть на правильном пути. Предполагая, что беты начать с:

betas <- c(1.2,1.3,1.6,1.4,2.2) 

Я хочу, чтобы свести к минимуму такую ​​функцию, что b0 <= b1 <= b2 <= b3 <= b4

f <- function(x) { 
    x1 <- x[1] 
    x2 <- x[2] 
    x3 <- x[3] 
    x4 <- x[4] 
    x5 <- x[5] 

loss <- (x1 - betas[1])^2 + 
     (x2 - betas[2])^2 + 
     (x3 - betas[3])^2 + 
     (x4 - betas[4])^2 + 
     (x5 - betas[5])^2  

    return(loss) 
} 

Чтобы показать, что функция работает, потеря должна быть равна нулю, если мы передаем оригинальные беты в :

> f(betas) 
[1] 0 

и относительно большой с некоторыми случайными входами:

> set.seed(42) 
> f(rnorm(5)) 
[1] 8.849329 

И минимизированы при значениях я был в состоянии вычислить в Excel:

> f(c(1.2,1.3,1.4,1.4,2.2)) 
[1] 0.04 

Answer 1

1. Поскольку цель квадратичный и ограничения линейны, вы можете использовать solve.QP.

Он находит b, что сводит к минимуму

(1/2) * t(b) %*% Dmat %*% b - t(dvec) %*% b 

при ограничениях

t(Amat) %*% b >= bvec. 

Здесь мы хотим b, что сводит к минимуму

sum((b-betas)^2) = sum(b^2) - 2 * sum(b*betas) + sum(beta^2) 
        = t(b) %*% t(b) - 2 * t(b) %*% betas + sum(beta^2). 

С последнего члена, sum(beta^2), постоянна , мы можем отказаться от него, и мы можем установить

Dmat = diag(n) 
dvec = betas. 

Ограничения являются

b[1] <= b[2] 
b[2] <= b[3] 
... 
b[n-1] <= b[n] 

т.е.

-b[1] + b[2]      >= 0 
     - b[2] + b[3]    >= 0 
       ... 
        - b[n-1] + b[n] >= 0 

t(Amat) так, что это

[ -1 1    ] 
[ -1 1    ] 
[  -1 1   ] 
[    ...  ] 
[    -1 1 ] 

bvec и равна нулю ,

Это приводит к следующему коду.

# Sample data 
betas <- c(1.2, 1.3, 1.6, 1.4, 2.2) 

# Optimization 
n <- length(betas) 
Dmat <- diag(n) 
dvec <- betas 
Amat <- matrix(0,nr=n,nc=n-1) 
Amat[cbind(1:(n-1), 1:(n-1))] <- -1 
Amat[cbind(2:n,  1:(n-1))] <- 1 
t(Amat) # Check that it looks as it should 
bvec <- rep(0,n-1) 
library(quadprog) 
r <- solve.QP(Dmat, dvec, Amat, bvec) 

# Check the result, graphically 
plot(betas) 
points(r$solution, pch=16) 

2. Вы можете использовать constrOptim таким же образом (целевая функция может быть произвольной, но ограничения должны быть линейными).

3. В целом, вы можете использовать optim если вы reparametrize проблемы в не-задачу оптимизации, например

b[1] = exp(x[1]) 
b[2] = b[1] + exp(x[2]) 
... 
b[n] = b[n-1] + exp(x[n-1]). 

Есть несколько примеров here или there ,

Answer 2

Хорошо, это начинает принимать форму, но все же есть некоторые ошибки. Основываясь на разговоре в чате с @Joran, кажется, я могу включить условие, которое установит функцию потерь на произвольно большое значение, если значения не в порядке. Это, похоже, работает, ЕСЛИ это расхождение происходит между первыми двумя коэффициентами, но не после этого. Мне сложно разобраться, почему это так.

Функция минимизации:

f <- function(x, x0) { 
    x1 <- x[1] 
    x2 <- x[2] 
    x3 <- x[3] 
    x4 <- x[4] 
    x5 <- x[5] 

loss <- (x1 - x0[1])^2 + 
     (x2 - x0[2])^2 + 
     (x3 - x0[3])^2 + 
     (x4 - x0[4])^2 + 
     (x5 - x0[5])^2  

    #Make sure the coefficients are in order 
    if any(diff(c(x1,x2,x3,x4,x5)) > 0) loss = 10000000 

    return(loss) 
} 

Работа пример (вид, кажется, потеря будет сведена к минимуму, если b0 = 1.24?):

> betas <- c(1.22, 1.24, 1.18, 1.12, 1.10) 
> optim(betas, f, x0 = betas)$par 
[1] 1.282 1.240 1.180 1.120 1.100 

Нерабочий пример (обратите внимание, что третий элемент по-прежнему больше второго:

> betas <- c(1.20, 1.15, 1.18, 1.12, 1.10) 
> optim(betas, f, x0 = betas)$par 
[1] 1.20 1.15 1.18 1.12 1.10