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 &lt; x &lt; 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.
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 &amp; x_0^{M-1} &amp; \cdots &amp; x_0 &amp; 1  \\
 x_1^M &amp; x_1^{M-1} &amp; \cdots &amp; x_1 &amp; 1  \\
 \vdots &amp; \vdots &amp; \vdots &amp; \vdots &amp; \vdots \\
 x_N^M &amp; x_N^{M-1} &amp; \cdots &amp; x_N &amp; 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.