# Mathematica Interpolation or approximation

I have lists of points X and Y.

{1212,1216,1217,1219...}, {1,2,3,4,5,6...}


I combine them to make a one table of points t = Thread[{x, y}]; Using that table I build graphic:

So I need to find function/formula out of this table of point. So the best variant as I think is using InterpolatingFunction[range, table] and InterpolatingPolynomial[data, var]. But since we have huge amount of data the polynomial will be really huge (if I'm correct for n points, polynomial will be n-1). So if I understand correct it is not possible to use interpolation for it. Right? Therefore, we need to use regression to find approximated function, using Fit[t, func, v];. To make a better result, as I see, better divide graph in half, where it shifted to the right (near x = 4000). And make two approximation graph instead of one.

Since I'm not using Mathematica in daily use (it is my first time in this program), I need help.

What function is best to use in Fit[t, fun, v]? And how can I build two approximation functions for one graph? Do I need to divide the initial graph by my self?

Are there other methods that I missed?

Here is what I tried to do:

x = Import["C:\...\x.txt","List"];
x = Sort[x];
l = Length[x];
y = Array[0 + # &, l];
funi = Interpolation[t];
InterpolatingFunction[{1, 200}, 0];
p = InterpolatingPolynomial[t, v]   //trying to get a function
{funi[1212]}
ListLinePlot[{t}]  //graphic


//----------------------------------- so since I don't see any way that interpolation will work, I decide work more on a way of approximation, but if you have a method that work with interpolation I would glad to see it. So here what I come to:

x = Import["C:\...\x.txt","List"];
x = Sort[x];
l = Length[x];
y = Array[1211 + # &, l];
l = Fit[t, {1, v, v^2}, v]
l1 = LinearModelFit[t, v1, v1];
Show[ListPlot[t], Plot[l1[v1], {v1, 0, 7000}]]


And I have

it gave me a function that I looking for:

848.726+0.362082*v-0.0000206136*v^2

however i would like to build approximation function, that go through as many points as it can. So how can I do that? As What I ideally trying to make is two approximation functions for each half of the graph that go through as many points as it can.

• Interpolation[data] will construction a piecewise interpolation (default order/degree 3). There are options to control the kind of interpolation and the order. Look it up and see if that's what you want. (Use Fit & its relatives, if you're trying to model the data and have some idea about what the model is or should be.) – Michael E2 May 18 '16 at 18:11
• @MichaelE2 Well yes you have option to control the intrpolationfunction, but you don't have this option when you use interpolatingPolynomial, and since I have more than 100 point mathematica give a really long polynomial with recursion depth more than 1024. Therefore, only Fit probably will work – alex_mike May 18 '16 at 19:09
• @alex_mike Don't use InterpolatingPolynomial, use the Interpolation command instead. Have you actually tried doing that? It should have no problem at all with a few thousand points. Fitting requires a model, and you clearly don't have one in mind. Interpolation is, I suspect, a much better option for you. – MarcoB May 18 '16 at 19:15
• @MarcoB Perhaps what alex_mike is trying to say is that he doesn't want a piecewise formula but a single formula. But really, it seems to me, one should first know whether it's appropriate to fit data or to interpolate it, however much data there is, and let that determine the tools to apply. – Michael E2 May 18 '16 at 19:21
• To help you decide, how about telling us where you got the data in the first place... – J. M. will be back soon May 18 '16 at 19:25

While using a computer often means you don't have to worry if there is a large number of polynomials approximating data piecewise, the OP wishes to find a simple polynomial or two that roughly approximates the data. Here is an approach. Please note that data fitting and smoothing is not my forte; but the mathematics used here is fun and too alluring for me not to want to share.

Since interpolation is considered, I'll assume the data represents a function. The goal, then, is to approximate this function. This approach, unlike Anton Antonov's, will rarely interpolate any of the data, but it will approximate it. I'll use Interpolation to represent the function. One can then approximate the function in whatever way, say by a series in orthogonal polynomials such as the Chebyshev polynomials. The advantage to using orthogonal polynomials is that the truncated series solves a certain least-squares approximation problem. The Chebyshev polynomials are convenient because the series is easy to compute.

Some data

img = Import["http://i.stack.imgur.com/AYJja.png"];
rgb = Binarize@*ColorNegate /@ ColorSeparate@img;
blue = Thinning@ImageSubtract[rgb[[1]], rgb[[3]]]


idata = Sort[Mean /@ SplitBy[PixelValuePositions[blue, 1], First]];
ListPlot@idata


Note: We took the mean of the second coordinates for a given first coordinate. Duplication of abscissae occurs in idata because it comes from an image. The actual data used by the OP to produce the image might not have this problem.

Partition the data

It will be hard to approximate data that exhibits broad oscillations (or changes in concavity). The OP seemed delete the horizontal segment in the middle. Let's do that. I suspect that breaking up the data into pieces will take manual intervention.

{{{x1}}, {{x2}}} = Position[idata, #] & /@
DeleteCases[SplitBy[idata, Last], x_ /; Length[x] < 8][[1, {1, -1}]]
(*  {{{148}}, {{188}}}  *)


Construct the series

We'll show the steps in the code below. The first step is to select which data set to process, the whole, the first part, or the last part.

(* Select data set (uncomment): *)
(*data=idata;*)           (* whole data set *)
data = idata[[;; x1]];    (* left segment *)
(*data=idata[[x2;;]];*)   (* right segment *)

{xdata, ydata} = Transpose@data;
xdata = data[[All, 1]];
xdom = MinMax[xdata];

prec = MachinePrecision;           (* precision *)
nn = 256;   (* max. degree: 16 should be sufficient; more here to illustrate convergence*)
cnodes = N[Cos[Pi*Range[0, nn]/nn], prec]; (* Chebyshev extreme points *)
yFN = Interpolation@data;
cc = Sqrt[2/nn] FourierDCT[        (* DCT-I solves for the Chebyshev coefficents *)
yFN /@ Rescale[cnodes, {-1, 1}, xdom], 1];
cc[[{1, -1}]] /= 2;                (* DCT-I is off by half at the end points *)


Now cc contains the Chebyshev series coefficients. Take any many as you please to form the Chebyshev series of degree n:

c[[ ;; n+1]] . ChebyshevT[Range[0, n], x]


What degree is good? (How many coefficients to use?)

It's a bit of an art to determine how to use the Chebyshev coefficients cc. The lower the degree the smoother the model; the higher the degree, the closer the approximation to the data. You can plot the polynomials for various degrees, calculate the sum of squares of the residuals at the data point and/or the so-called "energy": $$\int_a^b \left[f''(x)\right]^2\;dx$$ It's not real clear what the optimal solution would be. As the degree n increases, the sum of squares should converge to zero and the energy tend to increase. One might seek a balance between closeness of approximation (sum of squares) and smoothing (energy). Just how to weight each component is a matter of judgment (as far as I know). One might also just see how the sum of squares improves with each increase in degree.

Below I show as a measure of quality 1000 times the energy plus the sum of squares of the residuals. (Auxillary functions given below.)

Manipulate[
With[{d2cc = Nest[dCheb[#, xdom] &, cc[[;; n]], 2]},  (* 2nd derivative of Cheb. series *)
With[{cf = chebeval[cc[[;; n]]], dcf = chebeval[d2cc[[;; n - 2]]]},
Plot[cf@Rescale[x, xdom, {-1, 1}], {x, xdom[[1]], xdom[[2]]},
Prolog -> {Red, Point@data},
PlotLabel ->
Row[{"1000 Energy ", #1, " + SS ", #2, " = ", 1000 #1 + #2} &[
i0Cheb[timesCheb[d2cc, d2cc], {xdom[[1]], xdom[[2]]}],  (* integrate Cheb. series *)
N@Total[(ydata - cf /@ Rescale[N@xdata, xdom, {-1, 1}])^2]  (* sum of squares *)
]],
Frame -> True,
PlotPoints -> Min[Max[30, Ceiling[n/2]], Ceiling[(xdom[[2]] - xdom[[1]])/2]],
PlotRange -> {xdom, MinMax[ydata] + {-2, 2}},
]
]],
{{n, 10}, 3, nn, 1, Appearance -> "Labeled"}
]


One can get quite a close approximation to the data (SS < 0.17). Note this is different than the interpolating polynomial through these points. Even though the degree is 240, evaluation is numerically stable (chebeval uses Clenshaw's algorithm; see below).

Here is the entire data set with a degree-15 approximation:

Note: It would probably make sense to use the mean sum of squares, SS/Length[data] instead of SS, when comparing data sets. (I didn't bother at first because my initial focus was on an optimal approximant for a single data set.)

Auxiliary functions

The Manipulate code above uses the following to do some algebra, calculus and evaluation with Chebyshev series. These are all based on various recursive properties of the Chebyshev polynomials.

(* Multiplies two Chebyshev series *)
Clear[timesCheb];
timesCheb[aa_?VectorQ, bb_?VectorQ] := Module[{cc},
cc = ConstantArray[0, Length[aa] + Length[bb] - 1];
Do[
cc[[1 + Abs[i + j - 2]]] += aa[[i]] bb[[j]]/2;
cc[[1 + Abs[i - j]]] += aa[[i]] bb[[j]]/2,
{i, Length@aa}, {j, Length@bb}];
cc
];

(* Differentiate Chebyshev series *)
Clear[dCheb];
dCheb::usage =
"dCheb[c, {a,b}] differentiates the Chebyshev series c scaled over the interval {a,b}";
dCheb[c_] := dCheb[c, {-1, 1}];
dCheb[c_, {a_, b_}] := Module[{c1 = 0, c2 = 0, c3},
2/(b - a) MapAt[#/2 &,
Reverse@Table[c3 = c2; c2 = c1; c1 = 2 (n + 1)*c[[n + 2]] + c3,
{n, Length[c] - 2, 0, -1}],
1]
];

(* Complete definite integral of Chebyshev series *)
Clear[i0Cheb];
i0Cheb::usage = "i0Cheb[c, {a,b}] integrates the Chebyshev series over {a,b}";
i0Cheb[c_] := i0Cheb[c, {-1, 1}];
i0Cheb[c_, {a_, b_}] :=
Total[(b - a)/(1 - Range[0, Length@c - 1, 2]^2) c[[1 ;; ;; 2]]];

(* Returns function to evaluate a Chebyshev series (compiled version)
* Clenshaw's algorithm, c = {c0, c1,..., cn} Chebyshev coefficients
* The Dot product given above is just as good (c . ChebyshevT[Range[0,n], x])
*   when x is numeric.  *)
Clear[chebeval];
chebeval::usage =
"f = chebeval[c], y = f[x], evaluates Chebyshev series, x rescaled to -1 <= x <= 1";
chebeval[c0 : {__Real?MachineNumberQ}] :=
With[{c1 = Reverse[c0], deg = Length@c0 - 1}, Compile[{x},
Block[{c = c1, s = 0., s1 = 0., s2 = 0.},
Do[
s = c[[i]] + 2 x*s1 - s2;
s2 = s1;
s1 = s,
{i, deg}];
Last@c + x*s1 - s2], RuntimeAttributes -> {Listable}]];

• Thank you for your detailed and educational answer! When I worked on mine I did not see the requirement to have simple functions -- which is in the comments, not in the question itself. I tried using Chebyshev polynomials for my data with Quantile Regression, but did not get good results. – Anton Antonov May 20 '16 at 22:44
• @AntonAntonov Thanks! I enjoyed your answer and blog post on quantile regression. I hadn't heard of that before. Yes, the original question had some ambiguity. I hinted, perhaps too subtly, that I thought your approach is probably more appropriate. – Michael E2 May 21 '16 at 1:18
• @MichaelE2 guys this is really interesting answers... Thanks – alex_mike May 21 '16 at 1:55
• I don't remember chebfun being capable of this... ;) – J. M. will be back soon May 21 '16 at 13:05
• @MichaelE2 sorry for stupid question, Im just really just started, but how can I get resulted function? I mean something like in Anton Antonov post below. Your Auxiliary functions is working but I don't see result, and it is really hard for me on my level to understand everything. – alex_mike May 23 '16 at 4:27

[...] however i would like to build approximation function, that go through as many points as it can. So how can I do that? As What I ideally trying to make is two approximation functions for each half of the graph that go through as many points as it can.

This can be easily done with Quantile Regression.

### Data

First let us generate some data:

SeedRandom[124]
data = Abs@
LowpassFilter[
Accumulate@
Re@Fourier[
Table[RandomReal[{-.5, .5}] Sinh[
Exp[RandomReal[{-.5, .5}]^2]], {2^8}]], .4];
data = Transpose[{Range[Length[data]], data}];

data[[All, 1]] =
data[[All, 1]] +
RandomVariate[NormalDistribution[0, 0.07], Length[data]];

data = Join[data[[1 ;; 100]], data[[130 ;; -1]]];

ListPlot[data, PlotTheme -> "Detailed"]


Note that

• the x-axis points have random offsets from the regular grid, and

• the creation of two data parts by removing elements in middle.

### Clustering

Cluster the data into two parts:

dataCls = FindClusters[data, 2];
ListPlot[dataCls, PlotTheme -> "Detailed", PlotLegends -> Automatic]


### Quantile regression

Import["https://raw.githubusercontent.com/antononcube/MathematicaForPrediction/master/QuantileRegression.m"]


Fit quantile curves through the points (quantile 0.5):

qfuncs = QuantileRegression[#, 120, {0.5}, InterpolationOrder -> 2][[
1]] & /@ dataCls;


Note that we used a large number of B-splines for the regression fit.

Plot the approximation errors:

opts = {ImageSize -> Large, PlotRange -> All, PlotTheme -> "Detailed"};
Grid[{{ListLinePlot[
MapThread[{#1, qfuncs[[1]][#1] - #2} &, Transpose[dataCls[[1]]]],
opts],
ListLinePlot[
MapThread[{#1, qfuncs[[2]][#1] - #2} &, Transpose[dataCls[[2]]]],
opts]}}]


We can see that we get a good approximation for both parts. Moreover, when the plotted differences are 0 this means that the regression quantiles (curves) have passed exactly through the corresponding data points.

### The approximation polynomials

The function QuantileRegression returns piecewise polynomials. Higher or lower order polynomials can be obtained with InterpolationOrder. Here is how (part of) the regression quantiles look like after simplification:

qfexpr1 = Simplify[qfuncs[[1]][x]];
qfexpr2 = Simplify[qfuncs[[2]][x]];

Grid[{{qfexpr1, qfexpr2}}, Alignment -> Top]


### Using x-coordinates as knots

We can use the x-coordinates of the data as knots for the regression quantiles with linear interpolation:

qfuncs = QuantileRegression[#, #[[All, 1]], {0.5},
InterpolationOrder -> 1][[1]] & /@ detaCls;


Then we get these errors:

• , thanks for your answer, It is very helpful. However, in result I have huge amount approximation polynomials, so is it possibly to come to one or two polynomials which will describe approximation for the graph? I know that it is mean that it will increase possibility of error. – alex_mike May 19 '16 at 0:15
• @alex_mike The approximation is piecewise, at each interval there is only one low order polynomial. – Anton Antonov May 19 '16 at 0:54
• then is it possible to create approximation using one polynomial? I mean i just need find simple polynomial that describe the whole cluster. – alex_mike May 19 '16 at 1:04
• as you can see above my function is not complex, it probably can be described as cubic of logarithmic function – alex_mike May 19 '16 at 1:06