Function Repository Resource:

MonotonicInterpolation

Source Notebook

Make a monotonic InterpolatingFunction with continuous second derivative for monotonic data

Contributed by: Ted Ersek

ResourceFunction["MonotonicInterpolation"][data]

returns an InterpolatingFunction that has a continuous second derivative and is monotonic wherever the data is monotonic.

Details and Options

The data should have the form {{x1,y1},{x2,y2},} where the xi and yi are real numbers.
The data does not need to be sorted.
The data can have any Precision, but the InterpolatingFunction returned can only give MachinePrecision results.
The InterpolatingFunction returned approximates the smoothest possible function that is monotonic except for extremum near data samples that are local minima or maxima among the data.
The following options may be given:
"KnotsBetweenSamples"2number of knots inserted between each pair of data samples.
MaxIterationsAutomaticmaximum number of iterations performed in QuadraticOptimization
ToleranceAutomaticthe tolerance to use for internal comparisons in QuadraticOptimization
When the data is monotonic, there exists a monotonic InterpolatingFunction with continuous second derivative provided at least two appropriately-spaced knots are inserted between each pair of data samples (section 8.3 of this paper).
The time ResourceFunction["MonotonicInterpolation"][data] needs to evaluate is proportional to (Length[data]+2(Length[data]-1)"KnotsBetweenSamples")3. When given a long List of data, it may need an unacceptable amount of time.

Examples

Basic Examples (6) 

The data below is monotonic, and we get a monotonic interpolation:

In[1]:=
data = {{0.0, 0.0}, {0.7, 0.02}, {2.1, 0.06}, {3.0, 8.28}, {3.8, 57.6}, {5.1, 62.2}, {5.9, 62.24}, {7.0, 62.26}};
fun = ResourceFunction["MonotonicInterpolation"][data]
Out[2]=

Plot the function:

In[3]:=
Plot[fun[x], {x, 0, 7}, Epilog -> Point@data]
Out[3]=

Next we see the function fun perfectly interpolates the data:

In[4]:=
fun /@ Part[data, All, 1]
Out[4]=
In[5]:=
Part[data, All, 2]
Out[5]=

The following suggests the derivative of fun is positive over the interval {0,7}:

In[6]:=
Min@Table[fun'[x], {x, 0, 7, 0.0003}]
Out[6]=

The derivative of fun is well behaved:

In[7]:=
Plot[fun'[x], {x, 0, 7}, PlotRange -> All]
Out[7]=

The second derivative of fun is continuous:

In[8]:=
Plot[fun''[x], {x, 0, 7}, PlotRange -> All]
Out[8]=

Scope (3) 

MonotonicInterpolation can handle exact data provided all values are real numbers:

In[9]:=
data = {{0, E^-4}, {1/3, E^-3}, {1/2, 1/(1 + E)}, {3/4, 1/E}, {1, 4 - E}, {3/2, 2}};
f = ResourceFunction["MonotonicInterpolation"][data];
Plot[f[x], {x, 0, 3/2}, Epilog -> Point@data, PlotRange -> All]
Out[11]=

The InterpolatingFunction above only gives MachinePrecision results:

In[12]:=
f[1]
Out[12]=
In[13]:=
Precision@f[1]
Out[13]=

f[1] evaluates to a MachinePrecision approximation of 4-:

In[14]:=
f[1] == 4 - E
Out[14]=

The next data has a local maximum at {2.4,9} and the InterpolatingFunction will have a local max someplace between x= 2 and x= 4. MonotonicInterpolation tries to find the smoothest InterpolatingFunction with no local min and no other local max; more precisely, MonotonicInterpolation tries to minimize with only one extremum:

In[15]:=
data = {{0, 0}, {0.4, 0.9}, {1, 1.35}, {1.5, 1.4}, {2, 2.3}, {2.4, 9}, {4, 3}, {5, 0.3}, {5.5, 0.3}, {6, 0.2}, {7, 0.0}};
g = ResourceFunction["MonotonicInterpolation"][data];
Plot[g[x], {x, 0, 7}, Epilog -> Point@data, PlotRange -> All]
Out[17]=

The following suggests g' is positive for all values in the interval {0,2.84} and never positive for values in the interval {2.841,7}:

In[18]:=
FindRoot[g'[x] == 0, {x, 2.8}]
Out[18]=
In[19]:=
Min@Table[g'[x], {x, 0, 2.84, 0.00003}]
Out[19]=
In[20]:=
Max@Table[g'[x], {x, 2.841, 7, 0.00003}]
Out[20]=

Two adjacent data samples above are {5,0.3}, {5.5,0.3}. Since these samples have equal y-values, g' is zero at all x in the interval {5,5.5}:

In[21]:=
Table[g'[x], {x, 5.0, 5.5, 0.02}]
Out[21]=

Using the same data, but unsorted, results in the same InterpolatingFunction:

In[22]:=
data = {{1.5, 1.4}, {0.4, 0.9}, {7, 0.0}, {1, 1.35}, {0, 0}, {2, 2.3}, {5.5, 0.3}, {2.4, 9}, {4, 3}, {6, 0.2}, {5, 0.3}};
h = ResourceFunction["MonotonicInterpolation"][data];
Plot[h[x], {x, 0, 7}, Epilog -> Point@data, PlotRange -> All]
Out[24]=

Options (11) 

KnotsBetweenSamples (3) 

More knots between each pair of samples will give smoother interpolations:

In[25]:=
data = {{0, 0.0}, {1, 0.02}, {2, 0.4}, {2.5, 11.3}, {4, 12.02}, {5, 12.04}};
ListLinePlot[data]
Out[26]=

The derivative of the InterpolatingFunction changes more slowly as the setting of the KnotsBetweenSamples option is increased:

In[27]:=
(* Evaluate this cell to get the example input *) CloudGet["https://www.wolframcloud.com/obj/8c4e5251-e3e0-4216-a7cf-5fdba1143a2d"]
Out[28]=

The integral of the second derivative squared decreases as the setting of "KnotsBetweenSamples" is increased:

In[29]:=
NIntegrate[f2''[x]^2, {x, 0, 5}]
Out[29]=
In[30]:=
NIntegrate[f3''[x]^2, {x, 0, 5}]
Out[30]=
In[31]:=
NIntegrate[f4''[x]^2, {x, 0, 5}]
Out[31]=
In[32]:=
NIntegrate[f6''[x]^2, {x, 0, 5}]
Out[32]=

MaxIterations (4) 

MonotonicInterpolation passes the setting of MaxIterations directly to QuadraticOptimization:

In[33]:=
SeedRandom[1234];
xs = Accumulate@RandomReal[{0.2, 0.75}, 40];
ys = Accumulate@RandomReal[{0.02, 0.2}, 40];
data = Transpose@{xs, ys};
ListLinePlot[data]
Out[37]=

Below a low setting for MaxIterations allows a solution to be found in less time than the default setting, but the quality of the solution may be sacrificed:

In[38]:=
AbsoluteTiming[
 f10 = ResourceFunction["MonotonicInterpolation"][data, MaxIterations -> 10]]
Out[38]=
In[39]:=
AbsoluteTiming[
 default = ResourceFunction["MonotonicInterpolation"][data]]
Out[39]=

Integrating of the second derivative squared shows that default is smoother than f10. This is explained by the message above:

In[40]:=
NIntegrate[f10''[x]^2, {x, Min@xs, Max@xs}]
Out[40]=
In[41]:=
NIntegrate[default''[x]^2, {x, Min@xs, Max@xs}]
Out[41]=

There is a slight difference between the InterpolatingFunction objects found above:

In[42]:=
Plot[f10[x] - default[x], {x, Min[xs], Max[xs]}]
Out[42]=

Tolerance (4) 

The default setting has the Tolerance determined automatically, but it is not able to find an InterpolatingFunction in this next example:

In[43]:=
data = {{0.0, 1}, {0.06, 1}, {0.09, 5}, {0.0925, 5.02}, {0.15, 7}, {0.18, 10}, {0.25, 10}};
f = ResourceFunction["MonotonicInterpolation"][data]
Out[44]=

Two monotonic interpolating functions of the data above are found using specific settings for Tolerance:

In[45]:=
g = ResourceFunction["MonotonicInterpolation"][data, Tolerance -> 10^-6]
Out[45]=
In[46]:=
h = ResourceFunction["MonotonicInterpolation"][data, Tolerance -> 10^-5]
Out[46]=

We see there are significant differences between the above interpolating functions:

In[47]:=
Plot[{g[x], h[x]}, {x, 0.05, 0.19}, Epilog -> Point@data, PlotLegends -> "Expressions"]
Out[47]=

Integrating of the second derivative squared shows that g is smoother than h:

In[48]:=
NIntegrate[g''[x]^2, {x, 0, 0.25}]
Out[48]=
In[49]:=
NIntegrate[h''[x]^2, {x, 0, 0.25}]
Out[49]=

Applications (4) 

Given data representing the fraction of a population of dogs that died at or before a certain age, we use Interpolation to create a continuous function describing the fraction of dogs that that would have died by any age less than 22:

In[50]:=
data = {{0, 0}, {1, 0.007}, {3, 0.02}, {6, 0.07}, {8, 0.14}, {12, 0.23}, {15, 0.75}, {18, 0.91}, {20, 0.96}, {22, 0.99}};
fractionDied = Interpolation[data, Method -> "Spline", InterpolationOrder -> 3];
Plot[fractionDied[years], {years, 0, 22}, PlotRange -> {0, 1}, Epilog -> Point@data]
Out[52]=

The next cell indicates that 13.25% of dogs would die at or before 10.25 years age, while the data provided says 14.0% of dogs died at or before 8 years of age. These are not compatible, so the above is not a very good InterpolatingFunction for this application:

In[53]:=
fractionDied[10.25]
Out[53]=

The previous example requires a monotonically increasing InterpolatingFunction. Using the data above, we could have met that requirement if we used Interpolation[data,Method"Hermite"], but that is not guaranteed to give a monotonic InterpolatingFunction for every monotonic set of data. Instead we can use linear interpolation as in the next cell, and that is guaranteed to give a monotonic InterpolatingFunction for every monotonic set of data:

In[54]:=
fractionDied = Interpolation[data, InterpolationOrder -> 1];
Plot[fractionDied[years], {years, 0, 22}, PlotRange -> {0, 1}, Epilog -> Point@data]
Out[55]=

However, this linear interpolation doesn’t have a continuous first derivative. It is reasonable to expect that the actual curve is smooth, and MonotonicInterpolation below tries to find the smoothest possible monotonically-increasing InterpolatingFunction that interpolates the monotonic data. The interpolation shown below is monotonic over the whole range plotted and it’s first two derivatives are continuous:

In[56]:=
fractionDied = ResourceFunction["MonotonicInterpolation"][data];
Plot[fractionDied[x], {x, 0, 22}, PlotRange -> {0, 1}, Epilog -> Point@data]
Out[57]=

Properties and Relations (8) 

Below, MonotonicInterpolation returns an expression with the built-in head InterpolatingFunction:

In[58]:=
data = {{0.0, 0.0}, {1.1, 0.1}, {2.2, 0.4}, {3.4, 6.0}, {4.2, 6.8}, {5.5, 6.9}, {6.5, 7.2}};
fun = ResourceFunction["MonotonicInterpolation"][data]
Out[59]=

We can use any of the following methods on an InterpolatingFunction:

In[60]:=
fun["Methods"]
Out[60]=

Next we get a lists of the x-values and y-values used to define the InterpolatingFunction fun:

In[61]:=
xValues = fun["Coordinates"]
Out[61]=
In[62]:=
yValues = fun["ValuesOnGrid"]
Out[62]=

The InterpolatingFunction fun also has specified derivatives at each of the x-values. Those derivatives can be retrieved from the internal representation of the InterpolatingFunction shown below:

In[63]:=
InputForm[fun]
Out[63]=

Get the derivative of the fun[x] at corresponding x-values:

In[64]:=
derivatives = Part[fun, 4, 3, 2 ;; -1 ;; 2]
Out[64]=

Get an expression used in the definition of the InterpolatingFunction fun:

In[65]:=
expression = Transpose@{Transpose@xValues, yValues, derivatives}
Out[65]=

Next we see how the above expression was used to define the InterpolatingFunction fun:

In[66]:=
fun === Interpolation[expression]
Out[66]=

The InterpolatingFunction fun is effectively a different cubic polynomial over several intervals. If we evaluate fun[x] using an x outside the domain of fun, it uses the cubic polynomial form the nearest end of the domain. The next cell performs surgery on the above expression to make an new InterpolatingFunction that extends the domain of the original function using linear extrapolation:

In[67]:=
{{x1}, y1, yp1} = Part[expression, 1];
{{xn}, yn, ypn} = Part[expression, -1];
expression = Join[{{{-10000}, y1 + yp1 (x1 - 10000), yp1}}, expression, {{{10000}, yn + ypn (10000 - xn), ypn}}];
linearExtrapolation = Interpolation[expression];
Plot[{fun[x], linearExtrapolation[x]}, {x, -3, 10}, Epilog -> Point@data, PlotLegends -> "Expressions"]
Out[71]=

Possible Issues (2) 

MonotonicInterpolation will often fail to find an InterpolatingFunction when some samples are very close together as in this example:

In[72]:=
data = N@{{1, 1}, {2, 1.8}, {3, 5}, {3.006, 5.0001}, {5, 7}, {6, 10}, {7, 10}};
ListLinePlot[data]
Out[73]=
In[74]:=
f = ResourceFunction["MonotonicInterpolation"][data]
Out[74]=

The number of unknowns that QuadraticOptimization needs to determine is Length[data]+2(Length[data]-1)*"KnotsBetweenSamples", and the time required for QuadraticOptimization to find a solution is proportional to (number of unknowns)3. As a result, the time MonotonicInterpolation requires increases dramatically as the setting for "KnotsBetweenSamples" increases. Also, when the setting for "KnotsBetweenSamples" is too large, the limitations of MachinePrecision arithmetic often prevent QuadraticOptimization from finding a solution. An example of that is given in the next cell that takes a few minutes to evaluate:

In[75]:=
data = {{0, 0.0}, {0.7, 0.02}, {2.1, 0.06}, {3, 20}, {3.8, 20}, {5.1, 62.2}, {5.9, 62.24}, {7, 62.26}};
ResourceFunction["MonotonicInterpolation"][data, "KnotsBetweenSamples" -> 45]
Out[76]=

Neat Examples (1) 

MonotonicInterpolation is efficient enough to use in Manipulate when the List of data is not too long:

In[77]:=
Manipulate[
 data = {{0, 0.0}, {0.7, 0.02}, {2.1, 0.06}, {3, y1}, {3.8, y2}, {5.1,
     62.2}, {5.9, 62.24}, {7, 62.26}};
 fun = Quiet@ResourceFunction["MonotonicInterpolation"][data];
 Plot[fun[x], {x, 0, 7}, Epilog -> Point@data],
 {{y1, 20}, 0.06, 20}, {{y2, 20}, 20, 62.2}]
Out[77]=

Publisher

Ted Ersek

Version History

  • 1.0.0 – 10 August 2020

Source Metadata

Related Resources

License Information