Ottimizzazione ristretta delle funzioni personalizzate in R

Question

Ho un modello combinato complesso per il quale posso definire una probabilità in una funzione e ho bisogno di ottimizzare i parametri. Il problema è che i parametri vanno in tutte le direzioni se non sono limitati. Quindi, devo implementare una restrizione dei parametri, e quello proposto dal professore è che la somma dei valori dei parametri quadrati dovrebbe essere uguale a 1

Ho giocato intorno con sia la funzione optim() e nlm(), ma Non posso davvero ottenere quello che voglio. La prima idea era di usare i parametri n-1 e calcolare l'ultimo dal resto, ma questo non funziona (come previsto).

Per illustrare, alcuni dati giocattolo e funzione riflette il problema nucleo di quello che voglio realizzare:

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) 

Ciò comporta ovviamente in:

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

Chiunque un'idea su come attuare restrizioni sui parametri nei processi di ottimizzazione?

PS: Sono consapevole del fatto che questo codice di esempio non ha alcun senso. È solo a scopo dimostrativo.

---

Modifica: Risolto! - Vedi la risposta di Mareks.

Answer 1

Penso che lo Ramnath answer non sia male, ma ha fatto qualche errore. La correzione alfa dovrebbe essere modificata.

Questa è migliorata versione:

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

ho ottenuto risultati simili alla versione illimitata.

---

[EDIT]

actualy questo non funziona correttamente se il punto di partenza è sbagliato. E.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 

Dà negate di vera stima perché mod <- glm.fit(m.mat,dd$Y) stime coefficiente negativo di X.

Penso che questa stima del glm non sia del tutto corretta. Penso che dovresti stimare l'intercettazione come media dei residui Y-X*alpha.

Qualcosa di simile:

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 

Answer 2

è necessario fornire maggiori dettagli sul vincolo. se hai a che fare con somma di quadrati uguali a uno, un modo elegante per risolverli usando optim è quello di lasciare che i parametri entrino in modo ottimale non vincolati, e riparameli all'interno della tua funzione di ottimizzazione.

per illustrare il mio punto, nell'esempio si è detto in precedenza, è possibile ottenere l'ottimizzazione esecuzione apportando le seguenti modifiche al codice:

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) 

questo dovrebbe funzionare per più di 3 variabili pure. fammi sapere se questo ha senso e se hai altre domande.

Answer 3

Potrebbe essere un po 'più complicato di quello che vuoi, e non ho il tempo di elaborare i dettagli al momento, ma penso che tu possa ancora farlo. Supponiamo che tenuti tutti i parametri tra 0 e 1 (si può fare questo con L-BFGS-B) e mappare l'Optim() Parametri p e il vostro parametri reali p' come segue:

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 qualcosa di un po' come quella.