Eingeschränkte Optimierung von benutzerdefinierten Funktionen in R

Question

Ich habe ein kompliziertes kombiniertes Modell, für das ich eine Wahrscheinlichkeit in einer Funktion definieren kann, und ich muss die Parameter optimieren. Problem ist, die Parameter gehen in alle Richtungen, wenn nicht eingeschränkt. Daher muss ich eine Einschränkung für die Parameter implementieren, und die vom Professor vorgeschlagene ist, dass die Summe der quadrierten Parameterwerte gleich 1 sein sollte.

Ich habe mit der optim() und nlm() Funktion gespielt, aber Ich kann nicht wirklich bekommen, was ich will. Die erste Idee war, n-1 Parameter zu verwenden und die letzte aus dem Rest zu berechnen, aber das funktioniert nicht (wie erwartet).

Zur Veranschaulichung einige Spielzeug Daten und Funktion das Kernproblem zu reflektieren, was ich erreichen möchte:

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) 

Dies führt offenbar in:

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

jemand eine Idee, wie Einschränkungen zu implementieren über Parameter in Optimierungsprozessen?

PS: Mir ist bewusst, dass dieser Beispielcode überhaupt keinen Sinn macht. Es ist nur zu Demonstrationszwecken.

---

Bearbeiten: Gelöst es! - Siehe Mareks Antwort.

Answer 1

Ich denke, dass Ramnath answer ist nicht schlecht, aber er macht einen Fehler. Die Alphakorrektur sollte geändert werden.

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

Ich habe ähnliche Ergebnisse zu Ihrer uneingeschränkten Version:

Dies ist die Version verbessert.

---

[EDIT]

Actualy wird dies nicht richtig funktionieren, wenn Startpunkt falsch ist. 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 

Es gibt negate der wahren Schätzung, da mod <- glm.fit(m.mat,dd$Y) Schätzungen negative Koeffizienten der X.

Ich denke, dass diese Glm Re-Schätzung ist nicht ganz richtig. Ich denke, Sie sollten Schnitt als Mittelwert der Residuen Y-X*alpha schätzen.

Etwas wie:

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

müssen Sie weitere Details zu Ihrer Abhängigkeit angeben. Wenn es sich um eine Summe von Quadraten handelt, die gleich Eins sind, besteht eine elegante Lösung darin, die Parameter optimierungsfrei zu lassen und sie in der Optimierungsfunktion neu zu parametrisieren.

meinen Punkt illustrieren, in dem Beispiel, das Sie oben angegeben haben, können Sie die Optimierung, indem sie die folgenden Änderungen an Ihren Code ausgeführt bekommen:

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) 

dies auch für mehr als 3 Variablen funktionieren würde. lassen Sie mich wissen, ob das Sinn macht und ob Sie weitere Fragen haben.

Answer 3

Es mag ein bisschen kniffliger sein, als Sie wollen, und ich habe nicht die Zeit, um die Details im Moment zu erarbeiten, aber ich denke, Sie können das immer noch tun. Angenommen, gebunden Sie alle Parameter zwischen 0 und 1 (Sie können dies tun, mit L-BFGS-B) und die Optim Karte() Parameter p und Ihre Parameter p‘wie folgt:

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

oder etwas ein bisschen wie Das.