Fitting by the right functions.

Ivan Cnop
icnop@vub.ac.be

Linear regression

Limited scope

Linear regression will be handled as a special case of polynomial approximations below. 
In fact there is no reason to handle (first order) linear regression separately if contemporary softwares handle polynomial approximations just as easily.
Moreover few functions or sets of data are handled well by linear regression.

values = Table[{x, Abs[x]}, {x, -3, 3,  .2}];

ListPlot[values]

Fit[values, {1, x}, x]

Chop[%]

Plot[%, {x, -3, 3}]

Show[{ListPlot[values], %}]

Linear regression is sometimes useful

The emphasis of linear regression in linear statistics is dated since technology allows for much finer (non-linear) statistical analysis and simulations.
Some situations where linear regression leads to new results will be treated in a separate text. 
An example is guessing the capacity in a logistic model from initial data (see later).

Polynomial approximations

Basis functions

A polynomial is a combination of powers of  x with real coefficients. 
We use the following shorthand for a list of powers of the variable   x  :

(x^Range)[0, 6]

Fitting a polynomial to a function

values = Table[{x, Abs[x]}, {x, -3, 3,  .2}];

ListPlot[values]

poly = Fit[values, (x^Range)[0, 6], x]

Plot[poly,  {x, -5, 5}, PlotRange -> {-3, 5}, 
      AspectRatio -> Automatic]

Show[{%, ListPlot[values]}]

Remark

Further outside the range of the data, polynomials may exhibit unwanted behavior

poly = Fit[values, (x^Range)[0, 6], x]

Plot[poly,  {x, -20, 20}, PlotRange -> All]

Show[{%, ListPlot[values]}]

In particular we will have to replace polynomials by other standard functions for (bounded) periodic data.

Fitting a polynomial to graphical inputs.

Select some points in an empty plot and save these

empty = 
        Plot[{x, 0}, {x, -3, 3}, PlotRange -> {-3, 5}, 
          AspectRatio -> Automatic];

values = {{-2.90971, 2.80376}, {-2.6919, 
            2.91266}, {-2.51039, 3.02157}, {-2.32888, 
            3.16678}, {-2.14737, 3.27568}, {-1.96586, 
            3.42089}, {-1.78435, 3.4935}, {-1.63914, 
            3.5661}, {-1.49393, 3.6024}, {-1.23982, 
            3.5298}, {-0.985701, 3.45719}, {-0.80419, 
            3.38459}, {-0.658982, 3.31199}, {-0.550075, 
            3.16678}, {-0.404867, 2.98527}, {-0.368565, 
            2.84006}, {-0.29596, 2.65855}, {-0.223356, 
            2.47704}, {-0.150752, 2.25922}, {-0.0781476, 
            2.07771}, {-0.0418456, 1.8599}, {0.0307586, 
            1.71469}, {0.103363, 1.56948}, {0.212269, 1.35167}, {0.321176, 
            1.13386}, {0.430082, 0.952347}, {0.57529, 0.734535}, {0.684197, 
            0.589326}, {0.829405, 0.444118}, {0.974614, 0.407815}, {1.26503, 
            0.48042}, {1.48284, 0.625628}, {1.66435, 0.734535}, {1.84586, 
            0.843441}, {1.99107, 0.879743}, {2.17258, 1.02495}, {2.31779, 
            1.13386}, {2.463, 1.24276}, {2.60821, 1.24276}, {2.78972, 
            1.31537}, {2.93493, 1.20646}};

ListPlot[values, PlotRange -> {-3, 5}, 
      AspectRatio -> Automatic]

poly = Fit[values, (x^Range)[0, 6], x]

Plot[poly,  {x, -5, 5}, PlotRange -> {-3, 5}, 
      AspectRatio -> Automatic]

Show[{%, ListPlot[values]}]

Experiment

Reduce the degree or increase the degree  n  of the base polynomials in Range[0,n]  to see what happens.

Fourier approximations

Fitting to a function

ListPlot[Table[{x, Abs[x]}, {x, -Pi, Pi, Pi/12}]]

values = Table[{x, Abs[x]}, {x, -Pi, Pi, Pi/24}];

ListPlot[values]

poly = Fit[values, (x^Range)[0, 6], x]

Please note that coefficients of odd powers of  x  are very small. This is because we are dealing with an even function.

Chop[%]

Plot[poly,  {x, -Pi, Pi}, PlotRange -> {-5, 5}, 
      AspectRatio -> Automatic]

Show[{%, ListPlot[values]}]

We use the following shorthand for a list of cosines . 
Since we are dealing with an even function  we expect we will not need sine functions.

Cos[Range[0, 7] t]

trig =  Fit[values, Cos[Range[0, 7] t], t]

Plot[trig, {t, -3 Pi, 3  Pi}]

Show[{%, ListPlot[values]}, AspectRatio -> Automatic]

Fitting to a periodic graphical input

We only have to draw input in the interval   [-π,+π]  if we know the data are periodic.

Plot[0,  {x, -Pi, Pi}, PlotRange -> {-1, 3}]

Save inputs by Copy-Paste and name these

{{-3.20275, 1.41396}, {-2.81207, 1.51328}, {-2.50714, 
          1.73763}, {-2.28798, 1.94913}, {-2.03388, 
          2.1618}, {-1.74484, 2.27397}, {-1.40656, 
          2.31136}, {-1.03335, 2.26112}, {-0.677607, 
          2.04962}, {-0.372684, 1.83695}, {-0.118581, 
          1.6126}, {-0.0677607, 1.56352}, {-0.0169402, 
          1.51328}, {0.0497616, 1.47589}, {0.0672312, 1.45135}, {0.270513, 
          1.26439}, {0.542085, 1.13936}, {0.608787, 1.1265}, {0.880359, 
          1.02718}, {1.35521, 0.939545}, {2.10005, 0.877614}, {2.35415, 
          0.976937}, {2.60826, 1.1265}, {2.81154, 1.23868}, {2.964, 
          1.31346}, {3.11646, 1.38825}};

values = %;

tab = 
      ListPlot[values, 
        PlotRange -> {{-π, π}, {-1, 3}}]

poly = Fit[values, (x^Range)[0, 5], x]

wide = 
      Plot[poly, {x, -3 Pi, 3  Pi}, 
        PlotRange -> {-1, 3}]

Which is not periodic.
To obtain periodic approximations, we use combinations of sines and cosines.
Therefore have to start with the union of these trigonometric base functions:

trig = 
      Fit[values, Union[Sin[Range[0, 3] t], Cos[Range[0, 3] t]], 
        t]

Plot[trig, {t, -Pi, Pi}, 
      PlotRange -> {{-π, π}, {0, 3}}]

Plot[trig, {t, -3 Pi, 3 Pi}, 
      PlotRange -> {{-3 π, 3 π}, {0, 3}}]

Show[%, tab, 
      PlotRange -> {{-3 π,  3 π}, {0, 3}}]

Outside the data interval the behavior is competely different. Recall

Show[wide, tab]

Other examples

If moreover we know the function is even or odd, we can restrict data to the interval[0,Pi]

max = 5

f[x_] := If[x < π/2, x, π - x]

evenValues  =  
        Table[{x, f[x]}, {x, 0, π, π/48}];

family = Table[Cos[n x], {n, 0, max}]

family // TableForm

Fit[evenValues, family, x]

Check the correspondenceo the right half of the interval:

Plot[{%, f[x]},  {x, -π, π}, 
      PlotRange -> {-1, 3}, 
      PlotStyle -> {{Hue[0]}, {Hue[0.6]}}]

The same holds for an odd function:

max = 5

f[x_] := If[x < π/2, x, π - x]

oddValues  =  
        Table[{x, f[x]}, {x, 0, π, π/48}];

family2 := Table[Sin[n x], {n, 0, max}]

family2 // TableForm

Fit[oddValues, family2, x]

Plot[{%, f[x]},  {x, -π, π}, 
      PlotRange -> {-2, 2}, 
      PlotStyle -> {{Hue[0]}, {Hue[0.6]}}]

We can easily check that the derivative of an even function is odd (and conversely)

max = 9

oddValues1 =  Table[{x, f'[x]}, {x, 0, Pi, Pi/48}];

Fit[oddValues1, family2, x]

Terms with exponents multiple of four are irrelevant:

Chop[%]

Check again the correspondence on the right half of the interval:

Plot[{%, f'[x]},  {x, -π, π}, 
      PlotStyle -> {{Hue[0]}, {Hue[0.6]}}]

Most calculus texts have classical examples similar to these.

Humps and Wavelets

Hump-like data

All previous fitting techniques go wrong if we have data which are localised around a given value of the x variable, and nearly zero far from this value. 
This is the case for many distributions in statistics ans many data sets in digital signal processing. 
It occurs whenever data lie in a stochastic way around some ideal value (called central value or mean in different situations).
We start with some data having a deviation from zero. If we know they are symmetric about the  y-axis, we can limit ourselves to the right half line. 
Feed such data by mouse in this empty graph and perform copy - save:

empty = 
        Plot[{x, 0}, {x, 0, 3}, PlotRange -> {-1, 3}, 
          AspectRatio -> Automatic];

values = {{0.0602322, 2.44638}, {0.139929, 2.40653}, {0.219626, 
            2.40653}, {0.339172, 2.36668}, {0.418868, 2.34675}, {0.51849, 
            2.30691}, {0.618111, 2.22721}, {0.697808, 2.18736}, {0.797429, 
            2.10766}, {0.837278, 2.02797}, {0.916974, 1.94827}, {0.956823, 
            1.86857}, {1.03652, 1.70918}, {1.09629, 1.56971}, {1.21584, 
            1.35054}, {1.31546, 1.13138}, {1.41508, 0.952057}, {1.47485, 
            0.832512}, {1.5147, 0.752815}, {1.53463, 0.653194}, {1.55455, 
            0.573497}, {1.5944, 0.473876}, {1.61432, 0.374254}, {1.65417, 
            0.294557}, {1.73387, 0.234785}, {1.79364, 0.155088}, {1.87334, 
            0.155088}, {1.95304, 0.155088}, {2.03273, 0.115239}, {2.11243, 
            0.095315}, {2.19213, 0.0753908}, {2.27182, 0.0753908}, {2.35152, 
            0.0753908}};

ListPlot[values, PlotRange -> {-1, 3}, 
      AspectRatio -> Automatic]

We  already suspected that polynomial approximation gives unsatisfactory results:

poly = Fit[values, (x^Range)[0, 6], x]

Plot[poly,  {x, -5, 5}, PlotRange -> {-3, 5}, 
      AspectRatio -> Automatic]

Show[{%, ListPlot[values]}]

The same holds for (even) trigonometric approximations:

trig = Fit[values, Cos[Range[0, 6] x], x]

Plot[trig,  {x, -10, 10}, PlotRange -> {-1, 3}, 
      AspectRatio -> Automatic]

Show[{%, ListPlot[values]}]

A choice of useful even basisfunctions

family1 = Sin[Range[6] x]/x

wave = Fit[values, family1, x]

Plot[wave,  {x, -10, 10}, PlotRange -> {-1, 3}, 
      AspectRatio -> Automatic]

Show[{%, ListPlot[values]}]

We know that the functions in family1 above decay as fast as   1/x. 
This means that the approximation will have this asymptotic behavior. 
Since our data  are going to zero faster than this - in fact they are supposed to be zero beyond the original frame - we may have a better result if we choose basis functions that tend to zero faster than those in family1:

Another choice of useful basisfunctions

Start with a bell-shaped function

Plot[(E^(-(x^2))), {x, 0, 2}]

and rescale x-axis by dividing the exponent by natural numbers:

family2  =  (E^(-(x^2) Range[6]))

wave = Fit[values, family2, x]

Plot[wave,  {x, -10, 10}, PlotRange -> {-1, 3}, 
      AspectRatio -> Automatic]

Show[{%, ListPlot[values]}]

This is quite satisfactory for the given data.

Changing the mean

Is achieved by translating the data over the x axis until the central value is zero.

Odd sets of data

If the data are odd, we replace the base sets by the following odd functions:

family1 = Cos[Range[6] x]/x

family2  =  x (E^(-(x^2) Range[6]))

Table[
      Plot[family2[[j]],  {x, -4, 4}, PlotRange -> {-1, 3}], {j, 
        6}]

empty = 
        Plot[{x, 0}, {x, 0, 3}, PlotRange -> {-1, 3}, 
          AspectRatio -> Automatic];

values = {{0.0403066, 0.573496}, {0.0403066, 0.673117}, {0.0403066, 
            0.772738}, {0.0602308, 0.852435}, {0.080155, 
            0.932132}, {0.100079, 1.03175}, {0.100079, 1.13137}, {0.100079, 
            1.21107}, {0.120003, 1.29077}, {0.120003, 1.39039}, {0.139928, 
            1.49001}, {0.159852, 1.58963}, {0.179776, 1.66933}, {0.219624, 
            1.74902}, {0.239549, 1.82872}, {0.279397, 1.92834}, {0.319245, 
            2.04789}, {0.33917, 2.12758}, {0.398942, 2.20728}, {0.458715, 
            2.28698}, {0.538412, 2.28698}, {0.618108, 2.28698}, {0.697805, 
            2.24713}, {0.717729, 2.16743}, {0.717729, 2.08774}, {0.737654, 
            2.00804}, {0.737654, 1.92834}, {0.757578, 1.82872}, {0.777502, 
            1.7291}, {0.797426, 1.6494}, {0.81735, 1.54978}, {0.837275, 
            1.47009}, {0.857199, 1.39039}, {0.897047, 1.31069}, {0.936896, 
            1.21107}, {0.976744, 1.13137}, {1.01659, 1.05168}, {1.05644, 
            0.97198}, {1.09629, 0.892283}, {1.13614, 0.812587}, {1.19591, 
            0.712966}, {1.23576, 0.633269}, {1.29553, 0.553572}, {1.37523, 
            0.473875}, {1.45492, 0.414103}, {1.53462, 0.35433}, {1.61432, 
            0.314482}, {1.69402, 0.274633}, {1.77371, 0.214861}, {1.85341, 
            0.175012}, {1.93311, 0.155088}, {2.0128, 0.155088}, {2.0925, 
            0.135164}, {2.1722, 0.0953156}};

ListPlot[values, PlotRange -> {-1, 3}, 
      AspectRatio -> Automatic]

wave1 = Fit[values, family1, x]

Plot[wave1,  {x, -10, 10}, PlotRange -> {-3, 3}, 
      AspectRatio -> Automatic]

Show[{%, ListPlot[values]}]

wave2 = Fit[values, family2, x]

Plot[wave2,  {x, -10, 10}, PlotRange -> {-3, 3}, 
      AspectRatio -> Automatic]

Show[{%, ListPlot[values]}]

Quite satisfactory indeed!