2010-10-22 26 views
5

Tengo un modelo combinado complicado para el que puedo definir una probabilidad en una función, y necesito optimizar los parámetros. El problema es que los parámetros van en todas las direcciones si no están restringidos. Por lo tanto, necesito implementar una restricción en los parámetros, y la propuesta por el profesor es que la suma de los valores de los parámetros cuadrados sea igual a 1.Optimización restringida de funciones personalizadas en R

He estado jugando con las funciones optim() y nlm(), pero Realmente no puedo obtener lo que quiero. La primera idea fue usar parámetros n-1 y calcular el último del resto, pero esto no funciona (como se esperaba).

Para ilustrar, algunos datos de juguete y la función que reflejan el problema central de lo que quiero lograr:

dd <- data.frame(
    X1=rnorm(100), 
    X2=rnorm(100), 
    X3=rnorm(100) 
) 
dd <- within(dd,Y <- 2+0.57*X1-0.57*X2+0.57*X3+rnorm(100,0,0.2)) 

myfunc2 <- function(alpha,dd){ 
    alpha <- c(alpha,sqrt(1-sum(alpha^2))) 
    X <- as.matrix(dd[,-4]) %*% alpha 
    m.mat <- model.matrix(~X) 
    mod <- glm.fit(m.mat,dd$Y) 
    Sq <- sum(resid(mod)^2) 
    return(Sq) 
} 

b <- c(1,0) 
optim(b,myfunc2,dd=dd) 

Esto se traduce obviamente en:

Error: (subscript) logical subscript too long 
In addition: Warning message: 
In sqrt(1 - sum(alpha^2)) : NaNs produced 

Alguien una idea sobre cómo implementar restricciones sobre los parámetros en los procesos de optimización?

PD: Soy consciente de que este código de ejemplo no tiene ningún sentido. Es solo para fines de demostración.


Editar: ¡Resolvió! - Ver respuesta de Mareks.

+0

¿Has probado 'constrOptim'? – James

+0

@James, lo miré hace algún tiempo, pero no pude encontrar una manera de traducir nuestra limitación de una manera factible. Lo veré de nuevo, gracias al puntero. Una de las cosas también es que -afaik- constrOptim es incluso más lento que optim, y ya tenemos problemas graves de rendimiento con el código. –

+0

¿Cuántos parámetros? –

Respuesta

2

Creo que Ramnath answer no está mal, pero comete un error. La corrección alfa debe ser modificada.

Esta versión se mejora:

myfunc2 <- function(alpha,dd){ 
    alpha <- alpha/sqrt(sum(alpha^2)) # here the modification ;) 
    X <- as.matrix(dd[,-4]) %*% alpha 
    m.mat <- model.matrix(~X) 
    mod <- glm.fit(m.mat,dd$Y) 
    Sq <- sum(resid(mod)^2) 
    return(Sq) 
} 

b = c(1,1,1) 
(x <- optim(b, myfunc2, dd=dd)$par) 
(final_par <- x/sqrt(sum(x^2))) 

me dieron resultados similares a la versión sin restricciones.


[EDIT]

podía comprender esto no funcionará correctamente si el punto de partida es erróneo. MI.g

x <- optim(-c(1,1,1), myfunc2, dd=dd)$par 
(final_par <- x/sqrt(sum(x^2))) 
# [1] -0.5925 0.5620 -0.5771 

Da cambiar signo de la verdadera estimación porque coeficiente negativo mod <- glm.fit(m.mat,dd$Y) estimaciones de X.

Creo que esta reestimación de glm no es del todo correcta. Creo que deberías estimar el intercepto como la media de los residuos Y-X*alpha.

Algo así como:

f_err_1 <- function(alpha,dd) { 
    alpha <- alpha/sqrt(sum(alpha^2)) 
    X <- as.matrix(dd[,-4]) %*% alpha 
    a0 <- mean(dd$Y-X) 
    Sq <- sum((dd$Y-a0-X)^2) 
    return(Sq) 
} 

x <- optim(c(1,1,1), f_err_1, dd=dd)$par;(final_par <- x/sqrt(sum(x^2))) 
# [1] 0.5924 -0.5620 0.5772 
x <- optim(-c(1,1,1), f_err_1, dd=dd)$par;(final_par <- x/sqrt(sum(x^2))) 
# [1] 0.5924 -0.5621 0.5772 
+0

Eso es una coincidencia. Acabo de llegar al mismo resultado al mismo tiempo. Thx para la confirmación –

+0

Estoy parado detrás de ti. Booo! ;) – Marek

+0

Supongo que estoy un poco sorprendido de que esto funcione, ya que es apropiado para n parámetros con solo n-1 grados de libertad. Tal vez sea solo una dificultad para las estadísticas, no para la optimización per se. @Joris, solo por curiosidad, ¿qué pasa con mi intento de solución a continuación que no se escala bien? –

1

necesita proporcionar más detalles sobre su restricción. si se trata de una suma de cuadrados igual a uno, una forma elegante de resolverlo usando optim es dejar los parámetros que ingresan optim sin restricciones, y repararlos en su función de optimización.

para ilustrar mi punto, en el ejemplo que usted ha indicado anteriormente, se puede obtener la optimización del funcionamiento al hacer los siguientes cambios en su código:

myfunc2 <- function(alpha,dd){ 
    alpha <- alpha^2/sum(alpha^2); 
    X <- as.matrix(dd[,-4]) %*% alpha 
    m.mat <- model.matrix(~X) 
    mod <- glm.fit(m.mat,dd$Y) 
    Sq <- sum(resid(mod)^2) 
    return(Sq) 
} 

b = c(1,1,1) 
optim(b,myfunc2,dd=dd); 
ans = b^2/sum(b^2) 

esto iba a funcionar durante más de 3 variables. avíseme si esto tiene sentido y si tiene preguntas adicionales.

+0

Resultó que mis primeras objeciones, basadas en la consideración teórica, fueron bastante desalentadoras. Marek corrigió tu código, pero me pusiste en el camino correcto. +1 para eso, un agradecimiento y una disculpa por mi comentario inicial. PD: Tuve que "editar" tu respuesta para poder emitir un voto. Agregué un espacio, así que nada cambió. –

+0

supuse que sus parámetros eran todos positivos. La modificación de marek permite ambos signos, que creo que funcionan mejor para su caso. – Ramnath

+0

No tiene nada que ver con los letreros. Tu normalización es incorrecta, 'alpha^2' no suma 1. Comprueba por ejemplo' alpha = c (0.5,0.5) ', después de tu transformación obtienes' c (0.5,0.5) 'cuyas sumas de square son' 0.5'. – Marek

0

Puede ser un poco más complicado de lo que usted quiere, y no tengo tiempo para resolver los detalles por el momento, pero creo que todavía puede hacer esto. Supongamos que consolidó la totalidad de parámetros entre 0 y 1 (se puede hacer esto con L-BFGS-B) y un mapa de la optim() Los parámetros P y su parámetros reales p' de la siguiente manera:

p_1' = p_1 
p_2' = sqrt(p_2*(1-p_1'^2)) 
p_3' = sqrt(p_3*(1-(p_1^2+p_2^2)) 
... 
p_n' = 1-sqrt(sum(p_i^2)) 

o algo un poco como ese.

+0

He estado jugando con esas cosas también, pero no se adaptan bien a problemas más grandes. –

Cuestiones relacionadas