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

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

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

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

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 я не могу понять, как установить это в R. Я уверен, что я просто немного плотный, но мог бы использовать некоторые указатели справа направление!

3/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

Ответ 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]

i.e.,

-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, если вы репараметрируете проблему в проблему без ограничений, например

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

Есть несколько примеров здесь или там.

Ответ 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