Квадратичное программирование Настройка ограничения

Я пытаюсь настроить простую модель OLS с ограничениями на коэффициенты в R. Код ниже работает. Однако это демонстрирует

y = c + a1x1 + a2x2 + a3x3 с ограничением a1+a2 = 1

Я хотел бы изменить это ограничение на: a1*a2 - a3 = 0

Спасибо за вашу помощь!

РАБОЧИЙ КОД:

'''

    set.seed(1000) 
    n <- 20
    x1 <- seq(100,length.out=n)+rnorm(n,0,2)
    x2 <- seq(50,length.out=n)+rnorm(n,0,2)
    x3 <- seq(10,length.out=n)+rnorm(n,0,2)
    constant <- 100
    ymat <- constant + .5*x1 + .5*x2 + .75*x3 + rnorm(n,0,4)
    xmat <- cbind(x1,x2,x3)

    X <- cbind(rep(1,n),xmat) # explicitly include vector for constant
    bh <- solve(t(X)%*%X)%*%t(X)%*%ymat

    XX <- solve(t(X)%*%X)
    cmat <- matrix(1,1,1) 
    Q <- matrix(c(0,1,1,0),ncol(X),1) # a1+a2=1 for y = c + a1x1 + a2x2 + a3x3
    bc <- bh-XX%*%Q%*%solve(t(Q)%*%XX%*%Q)%*%(t(Q)%*%bh-cmat)

   library(quadprog)
   d <- t(ymat) %*% X
   Rinv = solve(chol(t(X)%*%X)) 
   qp <- solve.QP(Dmat=Rinv, dvec=d, Amat=Q, bvec=cmat, meq=1, factorized=TRUE)
   qp

   cbind(bh,qp$unconstrained.solution)
   cbind(bc,qp$solution)

'''


r constraints quadratic-programming
person Katy    schedule 22.11.2019    source источник

Ответы (3)


arrow_upward
0
arrow_downward

Предполагая, что проблема состоит в том, чтобы минимизировать || ымат - Х б || ^2 при условии b[2] * b[3] == b[4] мы можем заменить b[4], что даст неограниченную nls проблему, показанную ниже. b ниже — это первые 3 элемента b, и мы можем получить b[4], умножив два последних элемента b ниже вместе. Пакеты не используются.

fm <- nls(ymat ~ X %*% c(b, b[2] * b[3]), start = list(b = 0:2))
fm

давая:

Nonlinear regression model
  model: ymat ~ X %*% c(b, b[2] * b[3])
   data: parent.frame()
     b1      b2      b3 
76.9718  0.6275  0.7598 
 residual sum-of-squares: 204

Number of iterations to convergence: 4 
Achieved convergence tolerance: 6.555e-06

Чтобы вычислить b4

prod(coef(fm)[-1])
## [1] 0.476805

Примечание

Аналогичным образом исходная задача (для минимизации той же цели, но с исходным ограничением) может быть сведена к задаче без ограничений и решена с использованием nls посредством подстановки:

nls(ymat ~ X %*% c(b[1], b[2], 1-b[2], b[3]), start = list(b = 0:2))

давая:

Nonlinear regression model
  model: ymat ~ X %*% c(b[1], b[2], 1 - b[2], b[3])
   data: parent.frame()
      b1       b2       b3 
105.3186   0.3931   0.7964 
 residual sum-of-squares: 222.3

Number of iterations to convergence: 1 
Achieved convergence tolerance: 4.838e-08

Можно было бы даже перепараметрировать, чтобы сделать эту первоначальную проблему решаемой с помощью lm

lm(ymat ~ I(X[, 2] - X[, 3]) + X[, 4] + offset(X[, 3]))

давать

Call:
  lm(formula = ymat ~ I(X[, 2] - X[, 3]) + X[, 4] + offset(X[, 3]))

Coefficients:
       (Intercept)  I(X[, 2] - X[, 3])              X[, 4]  
          105.3186              0.3931              0.7964
person G. Grothendieck    schedule 22.11.2019

arrow_upward
0
arrow_downward

Г. Гротендик - спасибо за ответ. К сожалению, это не сработало для меня.

Я решил разобраться с лагранжианом, который оказался для меня слишком сложным.

Потом понял,

a1*a2-a3 =0 
a1*a2 = a3
ln(a1*a2)= ln(a3)
ln(a1) + ln(a2) -ln(a3) = 0

Это оставляет мне дополнительное ограничение, которое я могу решить с помощью пакета quadprog.

person Katy    schedule 10.01.2020

arrow_upward
0
arrow_downward

Возможно, вы можете попробовать код ниже, используя fmincon()

library(pracma)
library(NlcOptim)
# define objective function
fn <- function(v) norm(ymat- as.vector( xmat %*% v),"2")
# the constraint a1*a2 - a3 = 0
heq1 = function(v) prod(v[1:2])-v[3] 
# solve a1, a2 and a3 
res <- fmincon(0:2,fn,heq = heq1)

такой, что

> res$par
[1]  1.9043754 -0.1781830 -0.3393272
person ThomasIsCoding    schedule 10.01.2020