• Sign Up
• Login
• API
• FAQ
• Tools
• Archive
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.

Top