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:我知道这个示例代码完全没有意义.它仅用于演示目的.

---

编辑:解决了! - 见Mareks的回答.

Answer 1

我认为拉姆纳特的回答还不错，但他犯了一些错误。应修改 alpha 校正。

这是改进版本：

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

我得到了与您的无限制版本类似的结果。

[编辑]

实际上，如果起点错误，这将无法正常工作。例如

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