# How to fit a function to data so that the fit is always greater than or equal to the data?

b = nst[n_] :=
Length[NestWhileList[If[EvenQ[#], #/2, 3 # + 1] &,
n, # > 1 &]];
nn = 500;
With[{stps = Array[nst, nn]},
Table[Max[Take[stps, n]], {n, nn}]
]


I'm working with the following list and I am trying to find a fit so that it's always greater than the data rather then the normal fitting method used in the FindFit function:

FindFit[b, x + y*Log[z], {x, y}, z]


I like the ability to change the fitting model in the FindFit function but I can't figure out how to set it for what I want. Help would be appreciated.

• I don't think the code you've posted does what you intended (I get b being Null). Can you check it? Jun 19, 2016 at 21:52
• Your right I forgot to mention i needed to do some finagling. I created the list with the first function and then reset it equal to b. the original b = didn't work. Jun 19, 2016 at 22:04
• You might try FindMinimum where the function you're minimizing is the distance from the curve to the data point when the curve is above the data point and some large value (say, 500) when the curve is below the observed data point (i.e., a penalty for the curve being below the data). But without any mention of how the data is generated in a probabilistic manner, this is just a manipulation of the data without the ability to make statistical inferences.
– JimB
Jun 19, 2016 at 23:02
• You may also be interested in the techniques shown here: mathematica.stackexchange.com/questions/114864/… Jun 19, 2016 at 23:10
• I have a hard time following but that just means more learning. also it shows its not as easy as I hoped. Jun 20, 2016 at 3:12

Create the list b as you have shown.

nst[n_] := Length[NestWhileList[If[EvenQ[#], #/2, 3 # + 1] &, n, # > 1 &]]

b = With[{stps = Array[nst, nn]}, Table[Max[Take[stps, n]], {n, nn}]];


It looks like

ListPlot[b, PlotStyle -> Blue]


It is apparent that we want to locate the first point in each group of horizontal points and use that in the constraint.

Those points can be located as follows:

data = Transpose@Join[{Range[500], b}];
(* {{{1, 1}, {2, 2}, {3, 8}, ..., {500, 144}} *)


data is a list of {index, b} pairs.

Next locate the positions where there is a jump.

pos = Position[Differences[b], x_ /; x > 0] + 1


Build a list of constraints forcing the desired function to exceed the y value at those positions.

constraints =
Map[x + y*Log[#[[1]]] >= #[[2]] &, Extract[data, pos]]
(* {x + y Log[2] >= 2, x + y Log[3] >= 8, x + y Log[6] >= 9,
x + y Log[7] >= 17, x + y Log[9] >= 20, x + y Log[18] >= 21,
x + y Log[25] >= 24, x + y Log[27] >= 112, x + y Log[54] >= 113,
x + y Log[73] >= 116, x + y Log[97] >= 119, x + y Log[129] >= 122,
x + y Log[171] >= 125, x + y Log[231] >= 128, x + y Log[313] >= 131,
x + y Log[327] >= 144} *)


Use the constraints in FindFit.

solution =
FindFit[b, {x + y*Log[z], Sequence @@ constraints}, {x, y}, z]
(* {x -> 69.7139, y -> 12.8302} *)


Plot it to validate the solution

Show[
ListPlot[list, PlotStyle -> Blue],
Plot[Evaluate[x + y*Log[z] /. solution], {z, 1, 500},
PlotStyle -> Black]
]


• nice, it actually works fine if you brute force put every point into constraints Jun 22, 2016 at 21:42

I had the same idea as Jim Baldwin, as constraints are often implemented as penalty functions. Here is one that severely penalized any negative residual. The parameter scale might need to be adjusted to be a significant fraction of the range of the data values.

ClearAll[penalty];
penalty[residuals_?VectorQ, scale_: 10] :=
scale*Length@residuals*(1 - Sign@Min[residuals]);


Here's how to use it with FindFit:

scale = Max[b] - Min[b];
fit = FindFit[b, model = x + y*Log[z], {x, y}, z,
NormFunction -> (Norm[#1] + penalty[#1, scale] & ),
Method -> "NMinimize"]


It will warn that the norm is nonlinear and therefore nonlinear methods will be used.

Plot[model /. fit, {z, 1, Length@b},
Epilog -> {Red, Point@Table[{i, b[[i]]}, {i, Length@b}]}]


Check that the fitted model stays above the data:

model - u /. fit /. {z -> Range@Length@b, u -> b} // Min
(* 6.40443*10^-8  *)