Ograniczona optymalizacja niestandardowych funkcji w R

Question

Mam skomplikowany połączony model, dla którego mogę zdefiniować prawdopodobieństwo w funkcji, i potrzebuję zoptymalizować parametry. Problem polega na tym, że parametry przechodzą wszystkie kierunki, jeśli nie są ograniczone. Dlatego muszę wprowadzić ograniczenie parametrów, a proponowane przez profesora jest to, że suma kwadratów wartości parametrów powinna wynosić 1.

Bawiłem się z funkcją optim() i nlm(), ale Naprawdę nie mogę dostać tego, czego chcę. Pierwszym pomysłem było użycie parametrów n-1 i obliczenie ostatniego z pozostałych, ale to nie działa (zgodnie z oczekiwaniami).

Aby zilustrować niektóre dane zabawki i funkcja odzwierciedlające problem rdzenia, co chcę osiągnąć:

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) 

Wynika to oczywiście w:

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

Ktoś pomysł na temat sposobu wdrożenia ograniczeń o parametrach w procesach optymalizacji?

PS: Mam świadomość, że ten przykładowy kod nie ma żadnego sensu. Jest to tylko do celów demonstracyjnych.

---

Edytuj: Rozwiązano! - Zobacz odpowiedź Mareks.

Answer 1

Myślę, że Ramnath answer nie jest zły, ale popełnia błąd. Korekta alfa powinna zostać zmodyfikowana.

Jest to ulepszona wersja:

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

mam podobne wyniki do swojej nieograniczonej wersji.

---

[EDIT]

właściwie to nie będzie działać prawidłowo, jeśli punkt początkowy jest źle. 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 

Daje negować prawdziwego szacunku bo mod <- glm.fit(m.mat,dd$Y) szacunki ujemny współczynnik X.

Myślę, że ta ponowna ocena szacunku nie jest całkiem poprawna. Myślę, że należy oszacować przecinek jako średnią reszt Y-X*alpha.

Coś jak:

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

Musisz podać więcej informacji o swoim ograniczeniu. jeśli masz do czynienia z sumą kwadratów równą jednemu, eleganckim sposobem na jej rozwiązanie za pomocą optymalizacji jest umożliwienie wprowadzania parametrów optymalnych bez ograniczeń i reparametryzacja ich w ramach funkcji optymalizacji.

zilustrować mój punkt widzenia, na przykład masz podaną powyżej, można uzyskać optymalizacja działa poprzez następujące zmiany w kodzie:

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) 

to będzie działać dłużej niż 3 zmienne, jak również. daj mi znać, jeśli to ma sens i jeśli masz dodatkowe pytania.

Answer 3

Może to być nieco trudniejsze, niż chcesz, a ja nie mam czasu na dokładne wyliczenie szczegółów, ale myślę, że nadal możesz to zrobić. Załóżmy, że związanie wszystkich parametrów pomiędzy 0 a 1 (można to zrobić z L-BFGS-B) i odwzorować Optim() Parametry P i twoje prawdziwe parametry p”następująco:

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

lub coś trochę jak że.