Optimización con restricciones

Question

Estoy trabajando con la salida de un modelo en el que hay estimaciones de parámetros que pueden no seguir las expectativas a priori. Me gustaría escribir una función que fuerce estas estimaciones de utilidad en línea con esas expectativas. Para hacer esto, la función debe minimizar la suma de la desviación cuadrada entre los valores iniciales y las nuevas estimaciones. Ya que tenemos expectativas a priori, la optimización debe estar sujeto a las siguientes limitaciones:

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

Por ejemplo, las estimaciones de los parámetros primas a continuación están flipflopped para B2 y B3. Las columnas Delta y Delta^2 muestran la desviación entre la estimación del parámetro original y el nuevo coeficiente. Estoy tratando de minimizar la columna Delta^2. He codificado esto en Excel y mostrado cómo Solver de Excel optimizaría este problema proporcionando el conjunto de restricciones:

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 

Después de leer a través ?optim y ?constrOptim, no soy capaz de asimilar cómo configurar esto en R. Estoy seguro de que estoy siendo un poco denso, pero podría usar algunos indicadores en la dirección correcta.

24/3/2012 - Recompensa adicional porque no soy lo suficientemente inteligente como para traducir la primera respuesta.

Aquí hay un código R que debería estar en el camino correcto. Suponiendo que las betas comienzan con:

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

que desean minimizar la siguiente función tal que 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) 
} 

Para demostrar que la función trabaja, la pérdida debe ser cero si pasamos las betas originales en :

> f(betas) 
[1] 0 

y relativamente grande, con algunas entradas al azar:

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

y reducirse al mínimo en los valores que fue capaz de calcular en Excel:

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

Answer 1

1. Dado que el objetivo es cuadrática y lineal limitaciones, puede utilizar solve.QP.

encuentra el b que minimiza

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

bajo las restricciones

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

Aquí, queremos b que minimiza

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). 

Desde el último término, sum(beta^2), es constante , podemos soltarlo, y podemos establecer

Dmat = diag(n) 
dvec = betas. 

Las limitaciones son

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

es decir,

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

modo que t(Amat) es

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

y bvec es cero .

Esto nos lleva al siguiente código.

# 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. Puede utilizar constrOptim de la misma manera (la función objetivo puede ser arbitraria, pero las restricciones tienen que ser lineal).

3. De manera más general, se puede utilizar optim si Reparametrize el problema en un problema de optimización no restringida, por ejemplo

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

Hay algunos ejemplos here o there .

Answer 2

bien, esto está empezando a tomar forma, pero todavía tiene algunos errores. Basado en la conversación en el chat con @Joran, parece que puedo incluir un condicional que establecerá la función de pérdida en un valor arbitrariamente grande si los valores no están en orden. Esto parece funcionar SI la discrepancia ocurre entre los dos primeros coeficientes, pero no después. Tengo dificultades para analizar por qué ese sería el caso.

función a minimizar:

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) 
} 

Ejemplo de trabajo (una especie de, parece que la pérdida se minimiza si 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 

-no laborables ejemplo (tenga en cuenta que el tercer elemento es aún más grande que el segundo:

> 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