Estimate Taylor series coefficients from samples of a function

Say I have measurements of a function $y = y(x)$, sampled at $x_i$ with some noise, that could be approximated by a Taylor series expansion. Is there an accepted way of estimating the coefficients for that expansion from my measurements?

I could fit the data to a polynomial, but that's not quite right, because for a Taylor series the approximation should get better the closer you are to a central point, say x = 0. Just fitting a polynomial treats every point equally.

I could also estimate the various orders of derivatives at my point of expansion, but then I need to make decisions about what differentiating filters to use and how many filter coefficients for each. Would the filters for the different derivatives need to fit together somehow?

So does anybody know of established methods for this? Explanations or references to papers would be appreciated.

CLARIFICATION

In response to the comment below, my sampling is a rectangular window from an infinite function, that isn't necessarily band-limited but does not have strong high-frequency components. To be more specific, I'm measuring the variance of an estimator (measuring displacement in a medical ultrasound signal) as a function of a parameter of the estimator (the level of deformation or strain of the underlying tissue). I have a theoretical Taylor series for the variance as a function of deformation, and would like to compare it to what I get from simulation.

A similar toy example might be: say you have a function like ln(x), sampled at intervals in x with some noise added. You don't know what function it really is and you want to estimate its Taylor series around x=5. So the function is smooth and slowly varying for a region around the point you're interested in (say 2 < x < 8), but is not necessarily nice outside of the region.

The answers have been helpful, and some kind of least-squares polynomial fit is probably the route to take. What would make an estimated Taylor series different from a normal polynomial fit, though, is that you should be able to shave off higher-order terms, and have the polynomial still approximate the original function, just within a smaller range about your initial point.

So maybe the approach would be to do a linear polynomial fit using only data close to the initial point, followed by a quadratic fit with a little more data, cubic using a little more than that, etc.

• Some questions (which may or may not be relevant): By sampled, do you mean the function is/was band-limited to below some Fs/2 frequency? Are your samples a rectangular window of an infinite function, repeating function, or the complete function? Jan 5 '12 at 22:26
• Like Dilip pointed out in his answer, using a Taylor series expansion requires that you have knowledge of the function's derivative at all sample points. I suppose you could utilize your theoretical expression for the derivatives of $y(x)$, but that somewhat diminishes the utility of using an independent simulation to confirm your theory. In order to best emulate the Taylor series behaivor with respect to higher-order terms, an approach like what you suggested, using differing orders of polynomial fits, might be useful. Jan 9 '12 at 18:17

Instead of exact polynomial fitting, you could use a least-squares fit, which will find the polynomial of specified order that minimizes the total squared error between the fit and the measured $(x_i,y_i)$ pairs. This can help to mitigate the effects of noise on the fit.

Given measurements $y_i$ of a function $y = f(x)$ at domain values $x_i$ ($i = 0, 1, \ldots , N$), choose a polynomial order $M \le N$ (if $M = N$, then you're down to exact polynomial fitting, as the $N$ points uniquely determine an $M$th order polynomial). Then, set up a system of equations that are linear in the desired polynomial coefficients $p_k$:

$$y_i = p_Mx_i^M + p_{M-1}x_i^{M-1} + \ldots + p_1x_i + p_0, \;\;\;\;\;i = 0, 1, \ldots , N$$

The least-squares problem can be solved by arranging the measurements into matrix-vector form:

$$\mathbf{A} = \left [ {\begin{array}{cc} x_0^M & x_0^{M-1} & \cdots & x_0 & 1 \\ x_1^M & x_1^{M-1} & \cdots & x_1 & 1 \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ x_N^M & x_N^{M-1} & \cdots & x_N & 1 \end{array} } \right ], \;\;\mathbf{y}= \left [ {\begin{array}{cc} y_0 \\ y_1 \\ \vdots \\ y_N \end{array} } \right ]$$

The least-squares solution generates the vector of polynomial coefficients $\left[p_M, p_{M-1}, \ldots , p_0\right]$ that minimizes the total squared error in the above linear system. The solution can be computed as:

$$\mathbf{\tilde p} = (\mathbf{A^TA})^{-1} \mathbf{A^Ty}$$

It is worth noting that the matrix $(\mathbf{A^TA})^{-1} \mathbf{A^T}$ is also known as the pseudoinverse of the matrix $\mathbf{A}$. You can then use the least-squares polynomial coefficient vector $\mathbf{\tilde p}$ to evaluate the polynomial at any other $x$ values that you wish.

• In the case of equispaced abscissas, this is not different from applying Savitzky-Golay smoothing on your data.
– user276
Jan 11 '12 at 23:33
• Plus 1 for a nice answer. LSE is very ubiquitous indeed. Aug 8 '13 at 19:36

Ignore noise for now.

Given $n+1$ points $(x_i,y_i)$ where the $x_i$ are distinct numbers, you can, as you say, fit a polynomial $f(x)$ of degree at most $n$ through these points. Lagrange interpolation, for example, is a standard method for this. But, it is believed that the points are actually on a curve $y=g(x)$ where $g(x)$ is not necessarily a polynomial (e.g., it might be $e^x$ or $(x+a)/(x+b)$ etc.) and you would like to find the Taylor series for this function $g(x)$. Well, developing the Taylor series for $g(x)$ in the vicinity of $x=0$, say, requires knowledge of the value of $g(0)$ as well as the values of the derivatives $g^{(k)}(x) = \frac{\mathrm d^kg(x)}{\mathrm dx^k}, k=1,2,\ldots$ at $x=0$, while all that is known is the values of $g(x)$ at $n+1$ points $x_i$. Even if $x_i = 0$ for some $i$ so that $g(0)$ is known, it is still necessary to estimate $g^{(k)}(0)$ for $k=1, 2, \ldots$

Estimating the value of the derivatives of a function $g(x)$ at $x=0$ from its values $g(x_i)$ at selected points is a well-studied problem in numerical analysis, and the formulas to be used are readily available. What is not described in detail, or more commonly, not mentioned at all in the vicinity of these formulas, is that these formulas are obtained by fitting a polynomial $h(x)=\sum_k h_kx^k$ to the known points and estimating $g^{(k)}(0)$ as $h^{(k)}(0) = k!h_k$. Put another way,

From $n+1$ points $(x_i,g(x_i))$ of $g(x)$, we can develop the Taylor series for $g(x)$ only up to the term of degree $n$, and the truncated Taylor series is just $h(x)$, the polynomial that was fitted to the $n+1$ points.

So, what does fitting a polynomial mean? The standard fit is Lagrange interpolation which works well when there is no noise, the points $x_i$ are evenly spaced, and $0$ is the median value of the $x_i$. If noise is present, a least-squares fit of a polynomial of degree $m < n$ (see the answer by JasonR for details) is often better, and if we want to emphasize accuracy in the vicinity of $x=0$, a weighted least-squares fit can be used. Weighting the error terms from points in the vicinity of $0$ more than error terms from far away forces the minimization algorithm to produce an even better fit near $0$ at the expense of poorer accuracy far away from $0$. Of course, one also has to defend the choice of weighting function against naysayers who prefer a different weighting (or no weighting).

Example: Given $3$ points $(-1,y_{-1}), (0,y_0), (1,y_1)$, the Lagrange interpolation formula gives \begin{align*} f(x) &= y_{-1}\frac{x(x-1)}{2} -y_0(x^2-1) + y_{1}\frac{x(x+1)}{2}\\ &= y_0 + \frac{y_1-y_{-1}}{2}x + \frac{y_1-2y_0 + y_{-1}}{2}x^2 \end{align*} where the coefficients of $x$ and $x^2$ are the "three-point" formulas for the first and second derivative as given in Table 25.2 of Abramowitz and Stegun's Handbook of Mathematical Functions, that is, the Lagrange interpolation formula is the truncated Taylor series for a function $g(x)$ such that $g(-1) = y_{-1}, g(0) = y_0, g(1)=y_1$.