7
$\begingroup$ 
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.

Improve this question
$\endgroup$
6
  • $\begingroup$ I don't think the code you've posted does what you intended (I get b being Null). Can you check it? $\endgroup$
    – mikado
    Jun 19, 2016 at 21:52
  • $\begingroup$ 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. $\endgroup$
    – shai horowitz
    Jun 19, 2016 at 22:04
  • 1
    $\begingroup$ 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. $\endgroup$
    – JimB
    Jun 19, 2016 at 23:02
  • 4
    $\begingroup$ You may also be interested in the techniques shown here: mathematica.stackexchange.com/questions/114864/… $\endgroup$
    – MarcoB
    Jun 19, 2016 at 23:10
  • $\begingroup$ I have a hard time following but that just means more learning. also it shows its not as easy as I hoped. $\endgroup$
    – shai horowitz
    Jun 20, 2016 at 3:12

2 Answers 2

Reset to default
7
$\begingroup$

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]

Mathematica graphics

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

Mathematica graphics

Improve this answer
$\endgroup$
1
  • $\begingroup$ nice, it actually works fine if you brute force put every point into constraints $\endgroup$
    – george2079
    Jun 22, 2016 at 21:42
6
$\begingroup$

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

Mathematica graphics

Check that the fitted model stays above the data:

model - u /. fit /. {z -> Range@Length@b, u -> b} // Min
(* 6.40443*10^-8  *)
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.