# fitting a monotonic function?

Let's say I have a few points and derivatives given for a function and I'd like to find the best fit through the data, which is monotonic. For example, if I take a polynomial ansatz, I get

f[x_] = Sum[a[i] x^i, {i, 0, 4}];
constraints = f[x1] == y1 && f'[x1] == 0 && f[x2] == 0 && f[x3] == y3 && f'[x3] == 0;
result = f[x] /. Solve[constraints, Table[a[i], {i, 0, 4}]][[1]] // FullSimplify;


which for some xi and yi sometimes leads to an interpolation that is not monotonic:

as can be seen by the fact that the minimum is not at the region boundary.

How can I produce a fit (with the simplest model possible) that would satisfy all constraints and simultaneously return a result that is monotonic?

EDIT:

More details on the data:

The data is always given in terms of points at the far left and far right boundary of a range (positive on the left and negative on the right), slopes at these points are zero, and the position at which the function crosses the x-axis is also known, so it looks something like:

Example values might be:

Xleft=0;
Xright=1;
Xmid=0.41;


As a list of points:

fvalues={{0,1180},{0.41,0},{1,-570}};
fprimes={{0,0},{1,0}};


for which fourth degree polynomial fit produces the non-monotonic function plotted above. Basically, I'm wondering what would be the simplest function template to fit 3 given points and 2 given slopes to, such that the function is monotonic?

• Do you have data with measurement error? Do you only need the fit to be monotonic in a specific range? I think you need to give more details. – JimB Sep 8 '19 at 4:40
• @JimB I've added an edit to the question with some details. – Kagaratsch Sep 8 '19 at 11:01

Here is one possible approach:

(* build the function to fit *)
n = 4;
pol[x_] = Sum[a[i] x^i, {i, 0, n}]
(* a[0] + x a[1] + x^2 a[2] + x^3 a[3] + x^4 a[4] *)

fit = NonlinearModelFit[
Join[ (* the data - y=0 are values, y=1 are derivatives*)
Prepend[0] /@ fvalues,
Prepend[1] /@ fprimes
],
{
(* the model, which is pol[x] for y=0, and pol'[x] for y=1*)
(1 - y) pol[x] + y pol'[x],
(* the monotonicity constraint *)
Reduce[ForAll[x, Xleft <= x <= Xright, pol'[x] <= 0]]
},
Table[a[i], {i, 0, n}],
{y, x},
Weights -> Join[ (* weights for values and derivatives *)
Table[10, Length@fvalues],
Table[1, Length@fprimes]
]
]

Plot[
fit[0, x],
{x, Xleft, Xright},
Epilog -> {Red, Point@fvalues}
]


The idea is to construct a function with an additional parameter (y) which controls whether the model returns pol[x] or pol'[x]. This allows us to specify both function values and derivative values for the fit (since the scales of the two values are off, we need to specify some appropriate Weights). The biggest issue of this approach is that Reduce takes a while to simplify the constraint, especially for higher order polynomials. Additionally, the constraint will produce complex numbers for some inputs (due to the roots involved), which causes NonlinearModelFit to complain.

• Wow, this is interesting! I did not know NonlinearModelFit could take <= constraints on whole ranges of variable values. – Kagaratsch Sep 8 '19 at 13:14

It would seem that you have a strong belief that only 3 points, two zero slopes, and monotonicity will faithfully represent all of the values between. If you don't really know the form of the underlying function, don't want to use more the 3 points, and are using polynomials to approximate the function, then I think you're way optimistic.

I would start with a monotonic function right from the start and not have to worry about determining what relationships the coefficients of a polynomial must have. One such class of monotonic functions is a cumulative distribution function. Below I've used a linear function of a Beta distribution:

fvalues = {{0, 1180}, {0.91, 0}, {1, -570}};
{{xLeft, yLeft}, {xMid, yMid}, {xRight, yRight}} = fvalues;


Given that there are only 3 points there are an infinite combination of α and β values that will provide a monotonic function passing through the 3 points with slopes of zero at each end. If the midpoint is less than half way, then we set α to some value and solve for β. Otherwise we set β to some value and solve for α.

If[xMid < (xLeft + xRight)/2, α = 3; β =., α =.; β = 3];
sol = NSolve[yRight +
(yLeft - yRight)*(1 - CDF[BetaDistribution[α, β], (xMid - xLeft)/(xRight - xLeft)]) == 0,
If[xMid < (xLeft + xRight)/2, β, α]][[1, 1]]
(* β -> 5.67779 *)

(* Show the resulting function *)
Show[ListPlot[fvalues, PlotStyle -> {Red, PointSize[0.02]}],
Plot[yRight +
(yLeft - yRight)*(1 - CDF[BetaDistribution[α, β] /. sol,
(x - xLeft)/(xRight - xLeft)]), {x, xLeft, xRight}, PlotRange -> All]]


• You are right, there are a lot of different monotonic functions around so one really can't get around modeling the problem a bit more. I've looked more into it and turns out the model predicts the function to be given by Erf with points and slopes appropriately fitted. – Kagaratsch Sep 9 '19 at 3:06
• But doesn't the slope only equal zero at $-\infty$ and $\infty$ when using Erf? You describe a bounded range of values. One can use the logit-normal distribution function (which is bounded by zero or one but can be modified to any finite range). – JimB Sep 9 '19 at 3:58
• Right again! Slopes equal zero was a simplification for the question here on MA stackexchange, since I did not want to complicate it more than necessary. – Kagaratsch Sep 9 '19 at 12:57