This is a quadratic programming problem coming from constrained linear regression

$A b + \varepsilon = y$ and $B^T b >0 $

then minimize $\varepsilon^T \varepsilon=(A b - y)^T (A b - y)$ which is equivalent to $min(-d^T b + 1/2 b^T D b)$ where $D=2A^TA$, and $2y^TA=d^T$

I would like to get confidence interval in a quadratic programming problem. I am using R. For the curious I provide the below code where I would like to get the confidence interval of the estimated coefficients. I have tried bootstrapping but I am having issues with non-positive definite matrix, I am hoping another approach might be possible. Maybe just shock the estimates and calculate $min(-d^T b + 1/2 b^T D b)$ which is the cost function for quadratic program.

Amat = matrix(c(1,1,1,1,1,1,1,3,6,9,12,15,1,9,36,81,144,225,1,27,216,729,1728,3375), nrow=6)

Dmat = matrix( c(12,92,992,12152,92,992,12152,158600,992,12152,158600,2150552,12152,158600,2150552,29910872), nrow=4)
dvec = matrix(c(12,96.1128386,1086.2395754,13835.6202146), ncol=1)
Bmat = matrix( c(0,1,3,7,0,1,5,19,0,1,7,37,0,1,9,61,0,1,11,91,0,1,13,127,0,1,15,169,0,1,17,217,0,1,19,271,0,1,21,331,0,1,23,397,0,1,25,469,0,1,27,547,0,1,29,631,0,1,31,721,0,1,33,817,0,1,35,919,0,1,37,1027,0,1,39,1141,0,1,41,1261), nrow=4)
bvec = matrix(c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), nrow=1)

bEst <- solve.QP(Dmat, dvec, Bmat, bvec)

[1]  0.879426809  0.026242400 -0.005319406  0.00035462
  • 5
    $\begingroup$ You have nominally stated a deterministic quadratic programming problem, which is specified by the inputs Dmat, dvec, Amat, bvec. This problem has a deterministic solution. Where does the randomness come in? You need to provide a context. $\endgroup$ Nov 20, 2015 at 18:27
  • $\begingroup$ Thanks, I added further information. Let me know if it is more clear now $\endgroup$
    – adam
    Nov 23, 2015 at 11:15
  • 1
    $\begingroup$ I believe the difficulty here is that for components of $B^\prime b$ that are close to zero, the CI needs to terminate at zero itself. Related to this--and even thornier--is the problem of determining a CI for the components that have been forced to zero by the restriction. No simple analytical procedure is going to work and developing the asymptotics looks difficult. Consider a bootstrapped CI or something generic like that, assuming you have a much larger dataset. $\endgroup$
    – whuber
    Dec 1, 2015 at 16:39
  • $\begingroup$ this is my dataset. 6 points, yo = c(0.9978865,1.0054508,0.8994358,0.820207,0.9638653,1.3131546), b=c(1,3,6,9,12,15) $\endgroup$
    – adam
    Dec 1, 2015 at 16:47
  • $\begingroup$ In general one thing you could do is to run Monte Carlo simulation on a sub sample; iff you don't have time series data. In any case you have too few data. $\endgroup$
    – Cerbero
    Sep 5, 2016 at 15:39

1 Answer 1

yo = c(0.9978865,1.0054508,0.8994358,0.820207,0.9638653,1.3131546)
e <- yo-Amat%*%bEst$solution
R2 <- 1-crossprod(e,e)/crossprod(yo-mean(yo),yo-mean(yo))

degreesOfFreedom <- (nRow - polyOrder - 1)
paramVar <- diag(solve(Dmat*0.5))*crossprod(e, e)/degreesOfFreedom
bEst$paramVar <- paramVar
bEst$tValue <- bEst$solution/sqrt(paramVar)
bEst$pValue <- (1- pt(abs(bEst$tValue), df=degreesOfFreedom ))*2

I propose to use a solution similar to linear regression such as $(A^TA)^{-1}\sigma$ where $\sigma$ is error standard deviation. page 18 of nonlinear regression notes has similar formula

  • 2
    $\begingroup$ This answer looks incorrect: a student $t$ distribution is not the right one to use when any of the constraints holds (or is close to holding). Could you explain why you think this is correct? $\endgroup$
    – whuber
    Oct 8, 2019 at 13:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.