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. Das Problem ist, dass die Parameter in alle Richtungen gehen, wenn sie nicht eingeschränkt werden. Daher muss ich eine Beschränkung für die Parameter einführen, und die vom Professor vorgeschlagene ist, dass die Summe der quadrierten Parameterwerte gleich 1 sein sollte.

Ich habe sowohl mit dem optim() et nlm() Funktion, aber ich kann nicht wirklich bekommen, was ich will. Die erste Idee war, n-1 Parameter zu verwenden und den letzten aus den restlichen zu berechnen, aber das funktioniert nicht (wie erwartet).

Zur Veranschaulichung einige Spielzeugdaten und -funktionen, die das Kernproblem dessen, was ich erreichen möchte, widerspiegeln:

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 natürlich zu :

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

Hat jemand eine Idee, wie man Einschränkungen für Parameter in Optimierungsverfahren implementieren kann?

PS: Ich bin mir der Tatsache bewusst, dass dieser Beispielcode überhaupt keinen Sinn macht. Er dient nur zu Demonstrationszwecken.

---

Bearbeiten : Gelöst! - Siehe Mareks Antwort.

Answer 1

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

Dies ist die verbesserte Version:

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 wie bei Ihrer uneingeschränkten Version.

---

[EDIT]

Wenn der Startpunkt falsch ist, funktioniert das nicht richtig. Z.B.

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 ergibt das Negativ der wahren Schätzung, weil mod <- glm.fit(m.mat,dd$Y) schätzt einen negativen Koeffizienten von X .

Ich glaube, dass diese glm-Neuschätzung nicht ganz korrekt ist. Ich denke, Sie sollten den Achsenabschnitt als Mittelwert der Residuen schätzen Y-X*alpha .

Etwa so:

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

Wenn es sich um eine Summe der Quadrate gleich eins handelt, ist eine elegante Lösung mit optim darin zu sehen, dass die Parameter ohne Einschränkung in optim eingegeben werden und innerhalb der Optimierungsfunktion neu parametrisiert werden.

Um meinen Standpunkt zu verdeutlichen, können Sie in dem von Ihnen genannten Beispiel die Optimierung durch folgende Änderungen an Ihrem Code zum Laufen bringen:

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 würde auch für mehr als 3 Variablen funktionieren. Lassen Sie mich wissen, ob dies sinnvoll ist und ob Sie weitere Fragen haben.

Answer 3

Es ist vielleicht etwas komplizierter, als Sie wollen, und ich habe im Moment nicht die Zeit, die Details auszuarbeiten, aber ich denke, Sie können es trotzdem tun. Nehmen wir an, Sie binden alle Parameter zwischen 0 und 1 (Sie können dies mit L-BFGS-B tun) und bilden die optim()-Parameter p und Ihre realen Parameter p' wie folgt ab:

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 in der Art.