polyval

Evaluating Polynomial and Standard Error of Prediction for Polynomial Regression

Standard Error of Prediction

In the following, we use $\mathbb{K}$ to represent either the real line $\mathbb{R}$ or the complex plane $\mathbb{C}$. Consider the following polynomial regression model:

where $\epsilon_i$ is a random variable, $x_i\in\mathbb{K}$ and $p_1,p_2,\ldots,p_{n+1}$ are given polynomial coefficients. Let $y_i$ be an observation of $Y_i$ and $\hat{y}_i=f(x_i)$ be its predicted value. The prediction error is the difference $\hat{y}_i-y_i=f(x_i)-y_i$. We define the residual vector as

where

Let $\mathrm{SE}_i$ be the standard error of the prediction $y_i$. Since $\mathrm{SE}_i$ is undefined when $m\leq n+1$, we assume that $n+1< m$. Let $a_i$ be the $i$th row of $A$. For each $i$,

To improve numerical stability, we obtain the inverse $(A'A)^{-1}$ by QR decomposition with column pivoting $A=QRP'$, where $Q\in\mathbb{K}^{m\times m}$ is orthogonal, $R\in\mathbb{K}^{m\times (n+1)}$ and $P$ is an $(n+1)\times (n+1)$ permutation matrix. Note that

where $R_1\in\mathbb{R}^{(n+1)\times (n+1)}$ is upper triangular. For convenience, let $R_P=R_1P'$. Then, the numerator of $\mathrm{SE}_i$ can be obtained by

where $e_i=\left(a_iR_P^{-1}\right)'$ is a column vector. Therefore, the standard error of prediction is

This formula is used by polyfit to computer delta.

Usage

y = polyval(p, x)

• It evaluates a given polynomial at each element of x. Coefficients of the polynomial are given by p.

• x should be a numeric array.

• p should be a vector containing coefficents of the polynomial.

  • p(1) corresponds to the highest order term.
  • p(end) corresponds to the zeroth order term.
  • p should be non-empty.

• The output has the same size as x.

y = polyval(p, x, [], mu)

• x and p are as described in y = polyval(p, x) above.
• mu is a vector of two numbers, where mu(0) and mu(1) are, respectively, the mean and standard deviation of x.
• The output y gives values of the polynomial p evaluated at (x - mu(0)) / mu(1).
• y has the same size as x.

[y, delta] = polyval(p, x, S)

• x and p are as described in y = polyval(p, x) above.

• y contains values of the polynomial evaluated at x.

• delta contains the standard errors of the predictions y.

  • More specifically, delta(i) is the standard error of the prediction $y(i)$ at $x(i)$.
  • delta is computed using the last formula shown before.
  • delta has the same size as x.

• S should be a structure obtained by polyfit with fields R, normr and df.

  • S.R is $R_P$ discussed before.
  • S.normr is $\|Ap-y\|_2$, the 2-norm of the residual vector, as discussed before.
  • S.df is the degree of freedom defined by $m-(n+1)$.
  • If S.df is not positive, all entries in delta will be Inf.

[y, delta] = polyval(p, x, S, mu)

• mu is a vector of two numbers, where mu(0) and mu(1) are, respectively, the mean and standard deviation of x.
• y and delta are computed using the transformed $x$-coordinates (x - mu(0)) / mu(1).
• See [y, delta] = polyval(p, x, S) above for the description of the outputs y and delta.

Example

• Plotting a polynomial.

Input

r = [2 3 4 5];
% Creating a polynomial with roots r
p = poly(r);
x = linspace(1.8,5.2);
% Evaluating the polynomial
y = polyval(p,x);
plot(x,y,'-o','MarkerSize',2)

Output

• Fitting noisy data using a 2nd order polynomial and showing 95% prediction intervals. The $y$-coordinates of the data is generated by $y=x^2+\epsilon$, where $\epsilon\sim N(0,0.2)$.

Input

clear
% x and y are data points.
x=linspace(0,2);
y=x.^2+normrnd(0,0.2,[1 length(x)]);
% Get polynomial of degree 2 that fits the data.
[p,S]=polyfit(x,y,2)
scatter(x,y)
hold('on')

% Plot the polynomial and standard errors.
xx=0:0.01:2;
[yy,delta]=polyval(p,xx,S);
plot(xx,yy)
plot(xx,yy+2*delta,'k--')
plot(xx,yy-2*delta,'k--')
hold('off')

Output