R中自定义函数的约束优化

我想是的 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)) )

我得到了和你的无限制版本相似的结果。

[编辑]

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

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

为了说明我的观点,在上面的示例中,可以通过对代码进行以下更改来运行优化:

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个以上的变量。让我知道这是否合理,如果你还有其他问题。

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

或者类似的事情。