# Estimating a distribution using binned data

I have data from which information has been irretrievably lost through binning and I am trying to find a plausible distribution from which it derives. There is an academic paper suggesting that the underlying data is Zipf or Pareto distributed, but I believe that the paper is wrong. Just because one can fit a line to the log-log transform of survival data does not mean that the data is really Zipf or Pareto distributed!

The data (which you can see here, table 2b) is in an $n\times 2$ table in which the first column is a discrete range defining a "bin" and the second relevant column is an integer count. It thus states that there are y1 firms that have between min1 and max1 employees, that there are y2 firms that have between min2 and max2 employees, etc. I thus know that there are y1+y2+ ...+ yn firms that have at least min1 employees, that there are y2+y3 + ... + yn firms that have at least min2 employees, etc.

I would, of course, love to be able to use EstimatedDistribution with various guesses as to functional form and then use DistributionFitTest to see if the guess was plausible. But I can not figure out how to get Mathematica to find the best parameters for a distribution where all I have is binned data. I've floundered about with SurvivalModelFit, WeightedData, all to little avail. All help appreciated.

• You could try posting this at Cross Validated first to identify an appropriate method to tackle this problem... – sebhofer Jan 4 '14 at 9:33

I believe this is more a research problem than a Mathematica issue because EstimatedDistribution uses maximum likelihood estimation, not curve fitting. However, note that this built-in also accepts a ParameterEstimator option that allows you to use other methods, like the method of moments, more akin to curve fitting.

Inspired by this, I dare to propose the approach below. Although it may have statistical issues as the method of moments and others have.

As you pointed out, from the data, one really has a CDF. I am taking Table 2(b), bottom value of the interval and the column "Firms" to generate the following data:

data = {{1, 0.782054700833024}, {5, 0.915541200207107}, {10,
0.953737841741506}, {20, 0.976968965269283}, {100,
0.996148921112848}, {500, 0.999340326618122}, {750,
0.999559869493397}, {1000, 0.999666866923585}, {1500,
0.999771918945952}, {2000, 0.99982747114708}, {2500,
0.999861407705746}, {5000, 0.999930577761625}, {10000,
0.999964874580999}}


The first item in each pair is the bottom of the bracket, the second value is the sampled CDF based on the data.

Now, if you want to fit this to a Pareto distribution, you may try

FindFit[data, CDF[ParetoDistribution[k, a]][x], {k, a}, x]


which gives you {k -> 0.10358, a -> 0.668371}.

Now, being able to use DistributionFitTest is problematic because, actually, we don't have real data. So, I believe, any fit measurements or scores you would get with that would be of little value if, somehow, one generates a random sample out the binned data.

My suggestion would be to compare the estimated CDF with the binned CDF using some sort of distance.

Anyway, hope this helps.

A little bit of history

To illustrate my points made above, let me show you a previous attempt generating a random sample.

You can think of binned data as a spliced distribution of certain uniform distributions. If you have data formatted like

data = {{{1, 4}, 3705275}, {{5, 9}, 1060250}, {{10, 19},
644842}, {{20, 99}, 532391}, {{100, 499}, 88586}, {{500, 749},
6094}, {{750, 999}, 2970}, {{1000, 1499}, 2916}, {{1500, 1999},
1542}, {{2000, 2499}, 942}, {{2500, 4999}, 1920}, {{5000, 9999},
952}, {{10000, 100000}, 975}}


then the following builds such a distribution

MakeMixture[data_] := SplicedDistribution[
data[[All, 2]],
Append[data[[All, 1, 1]], Last[data] // First // Last],
UniformDistribution[{#[], #[] + 1}] & /@ data[[All, 1]]];


As you can see, I closed the top bracket which was "more than 10000". If you don't do that, you will have to speculate a distribution for the tail anyway - issue one. One way to close it is to consider the total labor force in USA, minus the first brackets and uniformly distribute that over the "more than 10000" firms, namely 975. You get an average of 129.5K employees for those 975 firms. I found that unlikely and arbitrarily closed the interval at 100K. That's issue number two.

You can then generate random samples of such distribution using RandomVariate. But beware that the top bracket is very unlikely, around 0.016% chance of failing into that. So, if you want a good sample, you may want to get more than 10K items. Then you would use EstimatedDistribution:

rawData = RandomVariate[dist, 10000];
EstimatedDistribution[rawData, ParetoDistribution[k, a]]


In one attempt I got ParetoDistribution[1.00138, 0.583712]. Very different from the parameters above. Was it right? In reality, rawData is not real data; just a resampling of an already problematic distribution; another issue. You could use DistributionFitTest against such rawData but IMHO lacks statistical value because it is not based on real data. Issue 4.

If you consider iterating this process through Monte Carlo simulation (resampling many times), are you going to average the parameters? I though that what you would do depends on how the underlying distribution behaves wrt its parameters... enter bayesian statistics (Issue 5 and I stopped counting). And even in this case, DistributionFitTest is not going to give you any sensible score averaging, IMHO.

Those were my previous thoughts.

• Nice- I did the same resampling method, results were wacko IMO, MC is probably a way to get 'closer', but as you've said this kind of bootstrapping depends on making up some kind of distribution for the bins. Interesting problem! OP - I'd wager the raw data is available - have you contacted the agency? – ciao Jan 2 '14 at 22:31
• IMHO this should be tackled using Bayesian estimation: take an uninformative prior, add binned data and derive a posterior dist for model parameters. I believe this is worth a research paper. – carlosayam Jan 2 '14 at 22:57

Here is a stab at the discrete Zipf problem:

v = {3617764, 1044065, 633141, 526307, 90386, 6060, 3038, 3044, 1533,
904, 1934, 975, 981};
bins = Partition[ {1, 5, 10, 20, 100, 500, 750, 1000, 1500, 2000,
2500, 5000, 10000, 100000}, 2, 1];
tot = Total@v;
binz[s_] := binz[s] =
tot Total@
Table[PDF[ZipfDistribution[100000, s], i] // N,
{i, #[], #[] - 1}] & /@ bins;


Try to fit using two different error measures:

sumsquares[s_] := Total@((binz[s] - v)^2);
sumsquaresrel[s_] := Total@((binz[s]/v - 1)^2);
fitsrel = s /. Last@FindMinimum[sumsquaresrel[s], {s, 0, 1}]
fits = s /. Last@FindMinimum[sumsquares[s], {s, 0, 1}]


(* .86, .47 *)

The resulting fits are "ok" , not great..

Show[{
ListLogPlot[v,Filling -> Axis],
ListLogPlot[binz[fitsrel], Joined -> True],
ListLogPlot[binz[fits], Joined -> True],
ListLogPlot[binz[(fits + fitsrel)/2], Joined -> True,
PlotStyle -> Dashed]}, PlotRange -> All] At a glance I'd say you have a piecewise distribution, with firms<~100 following one rule and larger firms another..

The standard approach for determining maximum likelihood with binned data is to determine the probability of an observation being in each bin given the functional form of the distribution. These probabilities are either found by summing individual probabilities of the discrete elements within a bin or by integrating from the lower to upper limit of each bin in the continuous case. Just Google "maximum likelihood binned data".

In the case of fitting a Zipf (a Zeta distribution) (which has a minimum value of 1 so I've removed the case of zero employees) here is some Mathematica code to do so.

(* Frequency counts *)

f = {3705275, 1060250, 644842, 532391, 88586, 6094, 2970, 2916, 1542,
942, 1920, 952, 975};


(* Lower and upper limits of bins *)

lower = {1, 5, 10, 20, 100, 500, 750, 1000, 1500, 2000, 2500, 5000, 10000};
upper = {4, 9, 19, 99, 499, 749, 999, 1499, 1999, 2499, 4999, 9999, 100000};


(* Direct log likelihood function *)

logL[s_] := N[Sum[f[[i]]*Log[Sum[1/(z^s*Zeta[s]),
{z, lower[[i]], upper[[i]]}]], {i, 1, Length[f]}]]


(* Much faster using some of Mathematica's built-in functions *)

n = Sum[f[[i]], {i, 1, Length[f]}]
logL2[s_] := Sum[f[[i]]*Log[HurwitzZeta[s, lower[[i]]] - HurwitzZeta[s, 1 + upper[[i]]]],
{i, 1, Length[f]}] - n*Log[Zeta[s]]

(* Find maximum likelihood estimate of s *)

NMaximize[{logL2[s], 1 < s < 10}, s]

(* Plot the log likelihood function *)

Plot[logL2[s], {s, 1.1, 2}]


One would still, of course, examine the observed counts against the expected counts from the maximum likelihood fit and estimate the standard error of the estimator.

On the comment that maximum likelihood has "major issues with bias". Yes, it can. But not always and not even "most of the time". One can simulate data from the fitted distribution and check to see if bias is a major problem for your particular distribution. Try it, you'll like it.

As suggested, this is more of a research question than a mathematica question, but if you implement some of the estimation in Mathematica and post it, I would love to see it.

Take a look at the empirics in the survey http://pages.stern.nyu.edu/~xgabaix/papers/pl-ar.pdf
A few important things to look at in this literature:

1. The standard estimation methods such as MLE and OLS have major issues with bias: http://pages.stern.nyu.edu/~xgabaix/papers/RankMinusOneHalf.pdf
2. It is difficult to distinguish between a Lognormal and a Pareto-tail when estimating from right-truncated data. Here is something from mathematica: http://demonstrations.wolfram.com/PowerLawTailsInLogNormalData/ but also search for maximum entry tests for power law tails, etc.

But I fear that noone will be very surprised if you conclude that it isn't a Zipf's law due to all of these econometric issues. I would be very careful with Bayesian methods in this specific case.