# Trying to fit an unknown extreme distribution

First up, I'm quite new to Mathematica so any hints on better code would be greatly appreciated.

I have some histogram frequency data from an unknown distribution that I'm trying to fit. Here's the histogram and how I draw it:

Needs["Histograms"]

(* Data in the form "rank count"*)
(* I have some missing ranks that I want to set as having a count of 0 *)
fd = Table[0, {i, Max[raw[[All, 1]]]}];
Do[fd[[raw[[i, 1]]]] = raw[[i, 2]], {i, Length[raw]}]
Histogram[Log[fd], FrequencyData -> True]


Needs["Histograms"] seems to be deprecated but I couldn't find a nice way of plotting frequency data without it.

The ranks are actually counts themselves, this is a histogram of the frequencies at which I've observed X number of things. Does that make sense? I'm slightly concerned that I'm confusing myself here :)

Now I have the data and the log plot seems to show a nice continuous curve I thought I could find a line to fit it. I followed these instructions: FindFIt with BinCounts, for using FindFit over frequency data, so far so good. Let's try a power law distribution:

centers = MovingAverage[raw[[All, 2]], 2];
counts = raw[[All, 1]];
centered = Table[{centers[[i]], counts[[i]]}, {i, Length[centers]}];

xmin = 1;
model = ((a - 1)/xmin)*(x/xmin)^(-a)
pars = FindFit[centered, model, {a}, x]
nlm = NonlinearModelFit[centered, model, {a}, x];
nlm[{"BestFit", "ParameterTable"}]


I'm not particularly expecting a power law distribution to work but it demonstrates my method. Here's what I get:

(-1+a) x^-a
{a->1.16913}
| Estimate | Standard Error | t-Statistic   | P-Value
----------------------------------------------------------
a   | 1.16913  | 122731.        | 9.52597*10^-6 | 0.999992
`

So, my questions. Is what I'm doing correct? If I managed to find a model that describes my data will FindFit estimate the parameters for me? Can anyone help with what that distribution might be? I've tried (grasping at straws):

• Zipf's law
• Power law
• Gumbel
• Laplace
• Frechet

Some have been reported as having a very good fit (P-value ridiculously low) but the lines don't really match up and when I plug in some numbers I get really bad results.

I feel like I'm being hopelessly naive in trying to do this. I've had my head stuck in these numbers for so long I don't really know what's going on :)

-
Cheers @kguler, thanks for putting the image in :) –  WilliamMayor Aug 9 '12 at 11:53
my pleasure. Welcome to Mma.SE. –  kglr Aug 9 '12 at 11:56
Do you happen to have version 8? Fitting distributions is a built-in there. –  Sjoerd C. de Vries Aug 9 '12 at 13:14

Your question strikes me more as a probability and statistics problem than a Mathematica problem, but I'll offer this extended comment.

You could create a list of every distribution Mathematica provides and a function to evaluate all of them. The following can give you some guidance of the method and its pitfalls:

Probabilistics - Monte Carlo

I'd suggest a different approach.

The Corp of Engineers has a good presentation that goes toward understanding your data. You can download it at:

Choosing a Probability Distribution

and Stat Trek has a great site for learning about probability and statistics:

Stat Trek

I've got a good flow chart to help narrow down distribution choices that I'll search for and post later if I find it.

All of the above just serves as background. I recommend that you ask yourself questions about the data you have.

• Is it discrete?
• Is it continuous?
• What kinds of processes generate it?
• What characteristics do they have?
• What else is it like?
• Does it have limits?
• Do you have enough data to represent the process?

All of these kinds of questions go to developing an explanation for why a given (or even some custom) distribution ought to model your observations.

Asking these kinds of questions get you much further in your understanding of the problem. Rather than getting stuck in the numbers the numbers may begin to reveal real insights into what you've begun to study.

-

Mathematica 8 knows an incredible amount of statistical distributions. From the WRI web site:

It has a number of functions to estimate distribution parameters from experimental data:

-