I fit a cubic function (in mathematica)

$$ y(x) = a + b x + c x^2 + d x^3 $$

to my data and obtained a function.

I have the error in each coefficient ($\sigma_a$, $\sigma_b$, $\sigma_c$, $\sigma_d$). Now since I need the model to predict $y$ values for a given $x$ value, I also need a corresponding error for the $y$ value (coming from the errors in the coefficients).

The error in the $x$ value (from our measurements) is negligible. I can't figure out the formula I would use when plugging in an $x$ value to determine its corresponding error in $y$.

I've exhausted my googling capabilities and my textbooks in error analysis. Any help is greatly appreciated.

  • $\begingroup$ you will also need to know the pairwise covariances between each of the coefficient estimates to get a variance for the predicted value of $y$. Assuming the predictor $x$ is not symmetric around 0, the covariance between, $x$ and $x^{2}$, for example, will be non-zero. From there (conditioned on the value of $x$), you just need to use the formula for the variance of a sum. $\endgroup$ – Macro Jun 18 '11 at 1:51
  • 1
    $\begingroup$ Polynomial regression seems like the obvious choice here, given that you are willing to consider you $x$ measurement error negligible. A "quick and dirty" approach (the first that came to my mind) is to do ordinary least squares of $y$ on $(1,x,x^2,x^3)$. Then just use standard prediction formula from that. Perhaps I am missing something? $\endgroup$ – probabilityislogic Jun 18 '11 at 3:15

The question appears to ask for a predicted value and for a prediction interval about that value.

The predicted value is obtained by means of the formula using the estimated parameters and a specified value of $x$:

$$\hat{y}(x) = \hat{a} + \hat{b} x + \hat{c} x^2 + \hat{d} x^3.$$

In general, the uncertainty in $\hat{y}(x)$ comes from three sources:

  1. Uncertainty in the parameter estimates $\hat{a}$, ..., $\hat{d}$ due to assumed randomness of the original data. An explicit probability model makes sense of this. In most settings, the model is additive. That is, for any datum $(x_i,y_i)$ it assumes

    $$y_i = a + b x_i + c x_i^2 + d x_i^3 + \varepsilon_i$$

    where $\varepsilon_i$ is a random variable with zero expectation. Furthermore, absent any information to the contrary, it is often the case that the $\varepsilon_i$ are independent of each other and have similar distributions. These assumptions allow us to estimate that common distribution from the data.

  2. Uncertainty in the actual value that would be observed at a given value $x$. This uncertainty is now directly evident in the preceding model: it is the contribution of $\varepsilon$.

  3. Uncertainty in the true value of $x$. In the original question this is explicitly assumed to be inconsequential.

Therefore we would expect the prediction interval formula to depend on three estimates: (i) the predicted value, (ii) the uncertainty in the predicted value due to the parameter uncertainty in (1) above, and (iii) the estimated variance of $\varepsilon$ (from (2) above). The parameter uncertainty plays a relatively small role for values of $x$ close to the average of the data $x_i$, but as $x$ moves away from this average, the parameter uncertainty makes an ever greater contribution.

Formulas can be found by searching this site for "prediction interval." In one thread, @Rob Hyndman gives a general formula that is directly applicable to this question. It assumes the parameters are estimated using least squares. For a design matrix $X$, the formula takes the form

$$\hat{y} \pm k_\alpha \hat{\sigma} \sqrt{1 + \mathbf{X}^* (\mathbf{X}'\mathbf{X})^{-1} (\mathbf{X}^*)'}.$$

($\mathbf{X}^*$ is determined by $x$ and $\mathbf{X}'\mathbf{X}$ encodes the variances and covariances of the $x_i$.)

This contains parts we can explicitly match to (i), (ii), and (iii):

(i) The interval is centered around the predicted value $\hat{y}$.

(ii) The parameter uncertainty appears in the term $\hat{\sigma} \sqrt{\mathbf{X}^* (\mathbf{X}'\mathbf{X})^{-1} (\mathbf{X}^*)'}$. The stuff inside the square root measures a squared distance between $x$ and the average of the $x_i$.

(iii) The contribution of $\varepsilon$ appears as the "$1$" (multiplied by $\hat{\sigma}$) added to the parameter uncertainty. We recognize this as the usual Pythagorean formula for the variance of a sum of independent random variables: take the square root of the sum of their squares.

As usual, the coefficient $k_\alpha$ is chosen to achieve a desired level of confidence. It depends on that level, $1-\alpha$, and on assumptions about how the $\varepsilon_i$ are distributed.

A comment to the original question prompts me to remark that considerations of correlation among the variables $1$, $x$, $x^2$, and $x^3$ are not relevant. Such correlations do of course exist in general--it is a rare situation in which our variables are mutually orthogonal--but the correlations are already accounted for in the least squares regression machinery. It is wise to check that the correlations are not so great that they introduce numerical instability in the solutions. That is unlikely in this particular problem.

Another comment also calls for another follow-up remark to emphasize one point: the presence of the second source of uncertainty (the contribution of $\varepsilon$, which is usually material) shows that we need to do more than consider the variances and covariances of the estimates $\hat{a}$, ..., $\hat{d}$: we also need to incorporate the contribution of $\varepsilon$ in the estimated uncertainty of $\hat{y}(x)$.

  • $\begingroup$ Really appreciate how meticulous and well informed your response is. I will respond more once I am sure I fully look through and understand your response (Statistics are definitely one of my weaknesses) $\endgroup$ – Illy Jun 20 '11 at 17:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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