SHARE
TWEET

Untitled

a guest Jan 17th, 2019 117 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. library(quadprog)
  2. # Simulating data from a linear model
  3. set.seed(2019)
  4. n <- 100
  5. x1 <- rnorm(n)
  6. x2 <- rnorm(n)
  7. y <- 0.3 * x1 - 0.3 * x2 + rnorm(n, sd=0.5)
  8. X <- cbind(x1,x2)
  9. # estimating a linear model with no constraints
  10. r <- lm(y ~ -1 + x1 + x2)
  11.  
  12. ## adding equality constraints
  13.  
  14. # e.g. the coefficients should sum up to b
  15. A <- matrix(1,2,1) # constraint matrix
  16. b <- 0.2
  17. # using solve.QP
  18. eqR <- solve.QP(t(X)%*%X, t(X)%*%y, A, b, meq=1)
  19. # manual implementation
  20. LHS <- rbind(cbind(t(X)%*%X, A), cbind(t(A), 0))
  21. RHS <- c(t(y)%*%X, b)
  22. eqM <- solve(LHS, RHS)
  23. # compare outcomes
  24. cbind(eqR$solution, eqM[1:2])
  25.          [,1]       [,2]
  26. x1  0.4450044  0.4450044
  27. x2 -0.2450044 -0.2450044
  28. # check constraint
  29. t(A)%*%eqR$solution; t(A)%*%eqM[1:2]
  30.      [,1]
  31. [1,]  0.2
  32.      [,1]
  33. [1,]  0.2
  34.  
  35. ## adding inequality constraints
  36.  
  37. # e.g. the coefficients should be larger than b
  38. A <- diag(2)
  39. b <- c(0.452,-0.33)
  40. ineqR <- solve.QP(t(X)%*%X, t(X)%*%y, A, b)
  41. # check constraints
  42. t(A)%*%ineqR$solution >= b
  43.      [,1]
  44. [1,] TRUE
  45. [2,] TRUE
  46. ## manual
  47. # here I was not able to reproduce manually ineqR
  48. # meaning solely using solve()
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top