Ограниченная оптимизация пользовательских функций в R

Question

У меня сложная комбинированная модель, для которой я могу определить вероятность в функции, и мне нужно оптимизировать параметры. Проблема в том, что параметры идут по всем направлениям, если не ограничены. Следовательно, мне нужно реализовать ограничение на параметры, и один предложенный профессором в том, что сумма квадратов значений параметров должен быть равен 1.

Я играл вокруг с как функции optim() и nlm(), но Я не могу получить то, что хочу. Первая идея заключалась в использовании параметров n-1 и вычислении последнего из остальных, но это не работает (как и ожидалось).

Для иллюстрации, некоторые игрушечные данные и функции, отражающие основную проблему того, что я хочу добиться:

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) 

Это приводит, очевидно, в:

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

Кто-нибудь представление о том, как реализовать ограничения по параметрам в процессах оптимизации?

PS: Я знаю, что этот примерный код не имеет никакого смысла. Это просто для демонстрационных целей.

---

Редактировать: решила! - См. Ответ Марекса.

Answer 1

Я думаю, что Ramnath answer не плохо, но он совершил ошибку. Альфа-коррекция должна быть изменена.

Это улучшенная версия:

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

Я получил аналогичные результаты в вашей неограниченной версии.

---

[EDIT]

Actualy это не будет работать корректно, если точка начала неправильно. E.г

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

Это дает свести на нет истинной оценки, потому что mod <- glm.fit(m.mat,dd$Y) оценки отрицательный коэффициент X.

Я думаю, что это пересмотр оценки glm не совсем корректен. Я думаю, вы должны оценить перехват как среднее значение остатков Y-X*alpha.

Что-то вроде:

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

вам необходимо предоставить более подробную информацию о вашем ограничении. если вы имеете дело с суммой квадратов, равной единице, элегантный способ ее решения с помощью оптимизации - позволить параметрам, входящим в режим оптимизации, неограниченно, а также репараметризировать их в пределах вашей функции оптимизации.

, чтобы проиллюстрировать свою точку зрения, в данном примере вы изложили выше, вы можете получить оптимизацию работы путем внесения следующих изменений в код:

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) 

это будет работать более 3-х переменных, а также. сообщите мне, если это имеет смысл, и если у вас есть дополнительные вопросы.

Answer 3

Это может быть немного сложнее, чем вы хотите, и у меня нет времени для разработки деталей на данный момент, но я думаю, вы все еще можете это сделать. Предположим, что вы связаны все параметры между 0 и 1 (вы можете сделать это с помощью L-BFGS-B) и карта Optim() параметры р и ваш реальный параметры р»следующим образом:

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

или что-то немного, как что.