fitting a monotonic function?

Question

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?

Answer 1

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.

Answer 2

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]]