2011-03-25 10 views
12

Estoy tratando de usar http://rss.acs.unt.edu/Rdoc/library/stats/html/constrOptim.html en R para hacer la optimización en R con algunas restricciones lineales dadas pero no puedo averiguar cómo configurar el problema.optimización restringida en R

Por ejemplo, necesito maximizar $ f (x, y) = log (x) + \ frac {x^2} {y^2} $ sujeto a restricciones $ g_1 (x, y) = x + y < 1 $, $ g_2 (x, y) = x> 0 $ y $ g_3 (x, y) = y> 0 $. ¿Cómo hago esto en R? Este es solo un ejemplo hipotético. No se preocupe por su estructura, en su lugar, estoy interesado en saber cómo configurar esto en R.

gracias!

+1

@G y otros: ¿alguien puede explicar por qué no se respeta la publicación cruzada? ¿Sería aceptable mencionar que estás publicando en forma cruzada con un enlace? No tengo sentimientos fuertes de una manera u otra, pero es probable que exista cierta claridad sobre el tema. Si esto es algo que la comunidad R en general ha tratado anteriormente, creo que vincularlo a esas discusiones sería un buen comienzo. – Chase

+0

Si realiza una búsqueda en la pestaña Meta para "publicación cruzada", encontrará una variedad de opiniones, la mayoría de las cuales son relativamente aceptables para la publicación cruzada. (La publicación cruzada simultánea, sin embargo, parece molestar a la mayoría de las personas). Existe una fuerte ética anti-publicación cruzada en los grupos de ayuda y primo que se establece en la Guía de publicación R-Help. Me cuesta mucho ver que la Guía de publicación sea probatoria en SO. –

+0

@Chase Cuando las personas cruzan la publicación y no le hacen saber a nadie que es muy posible que alguien escriba una solución para un problema que ya está resuelto, lo que puede ser una pérdida de tiempo. Personalmente, no me importa si las personas cruzan publicaciones siempre que lo den a conocer (y no es excesivo). – Dason

Respuesta

16

Configuración de la función era trivial:

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

la creación de la matriz de restricciones era problemático debido a la falta de mucha documentación, y recurrieron a la experimentación. La página de ayuda dice "La región factible está definida por ui% *% theta - ci> = 0". Así que probé y esto parecía "trabajo":

> 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 

por lo que poner en una fila para cada restricción/límite:

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 

Para este problema existe una posible dificultad de que para todos los valores de x la función va a Inf como y -> 0. Obtengo un máximo alrededor de x = .95 ey = 0 incluso cuando presiono los valores iniciales hacia la "esquina", pero sospecho que esto es no el verdadero máximo que hubiera supuesto que estaba en la "esquina". EDIT: Siguiendo esta razoné que el gradiente podría proporcionar "dirección" adicional y añade una función de gradiente:

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

Esto hizo "dirigir" la optimización de un poco más cerca de la c (0,999 ..., 0) esquina, en lugar de alejarse de ella, como lo hizo para algunos valores iniciales. Me quedo un poco decepcionado de que el proceso parece "la cabeza por el acantilado" cuando los valores de partida son cerca del centro de la región factible:

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 

Nota: Hans Werner Borchers registró un mejor ejemplo en I-ayuda que logró obtener los valores de las esquinas al establecer la restricción ligeramente alejada del borde:

> 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 
+0

Esto es muy útil. El ejemplo que publiqué no era el ideal, pero puedo ver cómo configurar la función. ¿Qué es una región factible? – user236215

+0

Una región factible es el conjunto de puntos que satisface todas las restricciones. Un triángulo en tu especificación. –

Cuestiones relacionadas