4
$\begingroup$

I am trying to find a fit to the cumulative distribution of a set of points using FindFit or NMinimize.

In particular, I would like to find the parameters of the cdf of the Beta Distribution that would minimize the uniform distance to the cumulative distribution of the above mentionned points.

So in particular, I am avoiding the use of NonLinearModelFit for reasons outlined here What is the difference between FindFit and NonlinearModelFit

However I get errors. I would appreciate a lot if someone could have a look at my code and let me know what am I missing.. (I tried a bunch of solutions neither works..)

So, here is the empirical cumulative probability function (in my problem it is a function derived from a non-parameteric kernel) but let us cook it up as a piecewise function:

 piece[x_] := Piecewise[{{x^3, 1 >= x >= 0}, {1, x > 1}}, 0]

Now I generate my "model" whose parameters I want to determine

funcr[al_?NumericQ, be_?NumericQ, x_?NumericQ] := 
 CDF[BetaDistribution[al, be], x];

So I can define first a norm which I will be later minimizing with respect to parameters $al$ and $be$ of the beta distribution

 norm[al_, be_, x_] := Abs[funcr[al, be, x] - piece[x]];
 max[al_, be_] := ArgMax[{norm[al, be, x], 0 <= x <= 1}, x]

And then I am trying to minimize this max function over all $\alpha$ and $\beta$ of the beta distribution (as given here by al_ and be_).

 NMinimize[{max[al, be], al >= 0, be >= 0}, {al, be}]

And at this stage I get the error:

 *The function value   
 \Abs[-0.028069671808523725`+funcr[al,be,0.04573179628771076`]] is not \
 a number at {x} = {0.04573179628771076`}. >>*

I know I am very pedestrian with my code, but I would really appreciate all your suggestions from which I can learn how to be a bit more sophisticated and especially correct....!

$\endgroup$
3
$\begingroup$

You could do for example:

int[al_?NumericQ, be_?NumericQ] := NIntegrate[(funcr[al, be, x] - piece[x])^2, {x, 0, 1}]
nm = NMinimize[{int[al, be], al >= 1, be >= 1}, {al, be}]
Plot[{piece@x, funcr[al, be, x] /. nm[[2]]}, {x, 0, 1}, 
     PlotStyle -> {{Thickness[.01], Red}, {Dashed, Thickness[.01], Blue}}]

Mathematica graphics

$\endgroup$
  • $\begingroup$ This is a great way to do it, I realized now what was my problem. Also I did not think of minimizing the sum of squared differences (that implies convergence in the max of the absolute difference, that I was interested in). Thanks a lot for your answer! $\endgroup$ – Kass Nov 6 '15 at 20:11
  • $\begingroup$ I have just one more question for this, I post it as a separate answer below, as it is too long for a comment... $\endgroup$ – Kass Nov 7 '15 at 0:36
3
$\begingroup$

I am not quite following what you are trying to do with norm or max.

The procedure I followed was to make some data from your fake empirical cumulative probability function.

piece[x_] := Piecewise[{{x^3, 1 >= x >= 0}, {1, x > 1}}, 0];
data = Table[{x, piece[x]}, {x, 0, 1, 0.02}];

Copy and paste your "model"

funcr[al_?NumericQ, be_?NumericQ, x_?NumericQ] := CDF[BetaDistribution[al, be], x]

In order to get an idea of a reasonable starting values I plot the data and the model using Manipulate.

Manipulate[
 Show[
  ListPlot[data, PlotStyle -> Black],
  Plot[funcr[al, be, x], {x, 0, 1}, PlotStyle -> Red],
  PlotRange -> All
  ],
 {{al, 2.0}, 0.001, 10, Appearance -> "Open"},
 {{be, 1.0}, 0.001, 10, Appearance -> "Open"}
 ]

Mathematica graphics

I actually could get close to the numerical solution by hand but let's allow FindFit to locate it.

Now run FindFit. We constrain the parameters to be greater than zero and provide starting values.

FindFit[data, {funcr[al, be, x], al > 0, be > 0},
          {{al, 2.0}, {be, 1.0}}, x]

(* {al -> 2.99998, be -> 0.999995} *)

Very quickly it converges on al = 3 and be=1.

You can go back to the Manipulate and validate that the data matches the model when the parameters are set to 3 and 1.

If you want to use something other than a least squares fit (i.e., Norm of 2) you can create it using the NormFunction option of FindFit.

For example here is the L1 norm

FindFit[data, {funcr[al, be, x], al > 0, be > 0},
     {{al, 2.0}, {be, 1.0}}, x, NormFunction -> (Norm[#, 1] &)]

I got the same answers using the L1 and L2 norms.

$\endgroup$
  • $\begingroup$ Thank you Jack! It works fine, and super clear all, my only worry is that the found parameters are not the global minimizers (as find fit apparently does not control for global min/max). That is why I was trying to use NMinimize of the distance, with the hope to get something "more global". In fact I am wondering how can one insure that the found parameters are indeed global minimizers of the Norm? Is there anything known on that? $\endgroup$ – Kass Nov 6 '15 at 20:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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