0

我想运行一个受约束的线性回归，使得系数是非负的并且总和为 1。通常这可以通过二次规划来完成，但是，我的参数（约束）比观察值（p > n）更多。我怎样才能做到这一点？

此 quadprog 解决方案仅适用于 n > p 的问题：

library(quadprog);
N <- 20
P <- 40
X <- matrix(runif(N*P), ncol=P)
btrue <- c(1,1,rep(0,(P-2)))
Y <- X %*% btrue + rnorm(N, sd=0.2)
C <- cbind(rep(1,P), diag(P))
b <- c(1,rep(0,P))
solve.QP(Dmat = t(X)%*%X, dvec = t(Y)%*%X , Amat = C, bvec = b, meq = 1)
4

1 回答 1

0

您可以使用以下模型执行此操作：

min r'r
r = X'b - y
b >= 0
sum(b) = 1
r free

其中 b 是要估计的参数，r 是残差。b 和 r 都是决策变量。这个问题总是凸的（即 Q=I 总是正定的）即使 n < p。

然而，Quadprog 可能仍然不喜欢它（我认为它需要一个严格的 pos def D 矩阵，不像大多数 QP 求解器）。通过将 D 更改为：

D = [ 0.0001*I  0  ]
    [ 0         I  ]

R 代码可能如下所示：

#
#
#  b = [b]  (P)
#      [r]  (N)
#
#  D =  [ 0 0 ] 
#       [ 0 I ]
#
#  d = [0]
#      [0]
#
#  A' =  [ 1' 0' ]
#        [ I  0  ]
#
# b0 = [1]
#      [0]
#


# fudge left upper sub matrix
D  = rbind( cbind( 0.0001*diag(P), matrix(rep(0,P*N),nrow=P,ncol=N)),
            cbind( matrix(rep(0,N*P),nrow=N,ncol=P), diag(N) ) 
           )
d = rep(0, P+N)

A = rbind( cbind( matrix(rep(1,P),nrow=1), matrix(rep(0,N),nrow=1)),
           cbind(diag(P),matrix(rep(0,P*N),nrow=P,ncol=N)))

b0 = rbind( matrix(c(1),nrow=1,ncol=1), matrix(rep(0,P),nrow=P,ncol=1))


solve.QP(Dmat = D, dvec = d , Amat = t(A), bvec = b0, meq = 1)
于 2020-04-07T18:01:58.690 回答