R 中自定义函数的受限优化

Question

我有一个复杂的组合模型，我可以在函数中为其定义似然，我需要优化参数。问题是，如果不受限制，参数会去各个方向。因此，我需要对参数进行限制，教授提出的参数值的平方和应等于1。

我一直在玩 optim()和 nlm()功能，但我无法真正得到我想要的。第一个想法是使用 n-1 参数并从其余参数中计算最后一个，但这不起作用(如预期)。

为了说明，一些玩具数据和功能反射(reflect)了我想要实现的核心问题:

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

我认为Ramnath answer不坏，但他犯了一些错误。应该修改 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