Lognormal distribution from world bank quintiles PPP data

I am not a stats person.

The world bank has data giving PPP (personal purchasing parity, or something like that) for quintiles (actually first 10%, second 10%, 2nd, 3rd, 4th 20%, 9th 10% and 10th 10%) of a country's population.

I've been reading a good deal about the Gini index and also modeling income distribution curves for populations, and I would like to turn this quintiles data into a lognormal distribution representing the same properties (e.g. that if China has 1,300 USD PPP for the 1st 10%, that after modeling it, the weighted average PPP (integration under the curve?) for the poorest 10% of the lognormal distribution would come out to 1,300).

Any thoughts on how to do this, tactically? I am cannot think my way through this - but I am a reasonable programmer of simple scripts and would use python's scipy and numpy to fit curves. Given some help on how to proceed.

@mpiktas, you are right, they do not give maximum income. The quintiles/decile information must be averages. I do not think you can straightforwardly fit this data as if it is raw data, to a distribution, if that were the case I would not have asked the question!

@Michael, given that the quintile/decile data points are average values for the entire quintile/decile bracket, is it possible still to come up with a best fit? Is this what you meant by least squares fitting?

-
10th decile would be maximum income? Does world bank really give such data, i.e. the maximum income of individual in population of specified country? –  mpiktas Sep 19 '12 at 12:43
I merged your accounts, TBrown, so you should be able to edit and comment freely from both computers. Cheers! –  whuber Sep 20 '12 at 1:56
@whuber - thank you!! –  TBrown Sep 20 '12 at 2:05
@Michael - Let me give a little context here because I think this is really fascinating stuff worth looking at. In China, a generic form of a medicine was priced at about 30% of the originator price, and achieved a 20x market share. I am wondering if this is because as you drop price, at a certain point along the income distribution (when going from rich to poor) you get a sharp increase in population - that is, when you've left the "long tail" of the rich and hit the bulk of the population. I want to see if I can account for this with the income distribution curve –  TBrown Sep 20 '12 at 2:12
Purchasing Power Parity, for the record –  Peter Ellis Dec 19 '12 at 11:32

A lognormal distribution is determined by two parameters, the mean and the variance of the related normal distribution. If you have raw data you could fit a lognormal distribution by maximum likelihood. If not you can use a fit criterion such as least squares or minimum sum of absolute errors to fit the given percentiles (quantiles) to values of a lognormal fit for these percentiles.

-
Dear all thank you for your posts - I am confused because the quintile data is, I think, PPP averaged over the quintile. So the data points aren't really points, exactly. Are they? fitting to an average seems like it will give me a totally wrong fit. Particularly for the 10th 10% of a country, which has to have a long, long tail. I can show you the data I want to model but I don't know how to add graphs like were added below. –  TBrown Sep 19 '12 at 16:07
This is the data: en.wikipedia.org/wiki/User:Pristino/… So, for example, China has: (0.1,0.2,0.4,0.6,0.8,0.9,1.0) for population brackets and (1367, 2444,3762, 5725, 8493, 12190, 24418) for average PPP for the bracket. It seems to me that if you graph Income vs. % population in the bracket, it looks a bit like a lognormal, but the left tail and right tail are far too high to fit. maybe I am thinking of this entirely the wrong way. –  TBrown Sep 19 '12 at 16:16

Here is the example of the quick and dirty R code to illustrate what Michael suggested:

Define quantiles available:

q<-c(0.1,0.2,0.4,0.6,0.8,0.9)


Create artificial data and add some noise

data <-jitter(qlnorm(q))


Create function to minimise

fitfun <- function(p)sum(abs(data-qlnorm(q,p[1],p[2])))


Run the optimiser with the initial guess of parameters of log-normal distribution:

opt <- optim(c(0.1,1.1))


The parameters fitted:

Display the fit visually:

aa<-seq(0,0.95,by=0.01)
plot(aa,qlnorm(aa,opt$par[1],opt$par[2]),type="l")
points(q,data)


Note, I intentionally plotted only 95%-quantile, since the log-normal distribution is unbounded, i.e. the 100%-quantile is infinity.

Usual caveats apply, real life example might look much uglier than this one, i.e. fit might be much worse. Also try Singh-Maddala distribution instead of log-normal, it works better for income distributions.

-
This works, but note it is not what Michael originally suggested (it is not the MLE); it is his alternative (fallback) suggestion. Your fitting criterion gives too much weight to the extreme quantiles relative to the middle ones (which might actually be a good thing in the circumstances, especially if one wants to emphasize income disparities). It also contains a somewhat arbitrary element: someone reporting the same data but providing more or different quantiles would change the nature of the fit by emphasizing precisely those quantiles they reported; the MLE doesn't have that problem. –  whuber Sep 19 '12 at 14:25
@whuber good points. Maybe the quantile fit would best be done by taking quantiles from say 1% and going to 99% in equal increments (maybe 5% increments) so that most of the distribution is represented. The other issue is the strong belief in the lognormal distribution. Goodness of fit to lognormal should be done before relying on the resulting fitted lognormal. –  Michael Chernick Sep 19 '12 at 16:24
Thank you for those thoughts, @Michael, especially the need to do GoF testing. Let me be more precise about my first point: the (possible) issues with quantile fitting are that (a) quantiles are correlated and (b) their variances change. This is why we can expect quantile fitting and MLE to produce different fits. Typically, quantile-fitting methods will not weight all quantiles equally (as shown in the example here). Instead, they weight quantile $q$ in proportion to $q(1-q)$ to compensate for the higher expected variability in the tails. –  whuber Sep 19 '12 at 17:22
@whuber I see. So you have argued the pros and cons of both approaches. Assuming the lognormal is a good fit to the data which approach would you recommend to the OP? –  Michael Chernick Sep 19 '12 at 17:28
@whuber, what would be likelihood function when you only have quantiles of distribution? Note Michael said fit MLE for raw data, in this case raw data is missing. –  mpiktas Sep 20 '12 at 6:50

I'm giving another answer, since more details about data were given. From the initial question it seemed that some quantiles are observed but that is not the case. The data is calculated in the following form.

1. Calculate the total income of all population
2. Divide population into income groups
3. Calculate the total income of population in groups defined in previous step.
4. Report for each group the proportion of total income in the group relative to total income of all population.

Suppose the population income is distributed according to unknown distribution function $F$. For the data the following income groups are defined:

1. Population with income in a range of $[0,F^{-1}(0.1))$
2. Population with income in a range of $[F^{-1}(0.1),F^{-1}(0.2))$
3. Population with income in a range of $[F^{-1}(0.2),F^{-1}(0.4))$
4. Population with income in a range of $[F^{-1}(0.4),F^{-1}(0.6))$
5. Population with income in a range of $[F^{-1}(0.6),F^{-1}(0.8))$
6. Population with income in a range of $[F^{-1}(0.8),F^{-1}(0.9))$
7. Population with income in a range of $[F^{-1}(0.9),\infty)$

For each of this range the following proportion is reported:

$$\frac{n_r\int_{l_r}^{u_r}xdF(x)}{N\int_{0}^{\infty}xdF(x)},$$

where $n_r$ is the number of people in a range $[l_r,u_r)$ and $N$ is the total population. The nominator of the fraction is number of people in a range times the average income in the range. Denominator has total number of people in the range times the average income.

Since the ranges defined are quantiles, proportions $n_r/N$ are known, i.e. for the first two and the last two ranges the proportion is equal to 0.1, for the rest 0.2.

The integral in nominator can be expressed in more convenient form:

$$\int_{l_r}^{u_r}xdF(x)=\int_{F(l_r)}^{F(u_r)}F^{-1}(u)du$$

The most obvious way to fit the data would be to integrate $F^{-1}$ numerically to a given range (or calculate the integrals analytically, which might be a challenge). Then calculate the proportions and fit them using your criterion of choice, least squares, least absolute deviations, etc. Note that one proportion is redundant since the proportions sum to one. Another caveat is that you need to know average income of the population, which is not given in the data.

-
I was thinking this way too, but there's a subtle problem: what is reported is not based on the quantiles of $F$ (we don't know it, so how could we possibly know its quantiles?), but rather it's based on the quantiles of the empirical distribution. For large data sets we could argue there's not much difference and maybe that's the way to go. In general, though, we can sort the data $x_{1}\le x_{2}\le\cdots\le x_{n}$ and what is reported are values of the form $(x_{r+1}+x_{r+2}+\cdots+x_{s})/(s-r)$ for suitably chosen pairs of ranks $(r,s)$. –  whuber Sep 20 '12 at 13:36
If you elect to use the large data set approximation, this particular problem simplifies greatly, because the truncated expectations for the lognormal distribution can be expressed in terms of the normal CDF $\Phi$ and its inverse: numerical integration is unnecessary (which makes finding the MLE much more practicable). –  whuber Sep 20 '12 at 13:38
@whuber, well the data description says very clear that it is based on quantiles. But naturally these quantiles are quantiles of empirical distribution. We then make a hypothesis that the empirical distribution is approximated well by some distribution $F$. But what would be MLE for such problem? MLE is usually defined when we have sample. Here we have some observed functionals of the sample. –  mpiktas Sep 20 '12 at 13:54
There nevertheless is a likelihood for these data, that likelihood can explicitly be related to a lognormal model for the underlying values (although the relationship is messy), and so by maximizing this likelihood we can identify the parameters in the usual way. –  whuber Sep 20 '12 at 20:00

A log-normal distribution is fully defined by the pair of parameters $\mu$ and $\sigma$. Since you want to fit this distribution to your data, it's sufficient to estimate these two values. Normally, you would have access to the raw data, and would apply the standard the maximum likelihood estimators (MLEs) for $\mu$ and $\sigma$, which are straightforward: $$\hat{\mu} = \frac{1}{n}\sum_i \ln(y_i) = \langle \ln y \rangle\\ \hat{\sigma}^2 = \frac{1}{n}\sum_i (\ln(y_i)-\hat{\mu})^{2} \enspace .$$ That is, $\mu$ is the mean of the logarithm of your observed data $\{y_i\}$, and $\sigma$ is the standard deviation of the logarithm of the data.

But in this case, you don't have the raw data. Instead, you have some sketchy information about the cumulative distribution function (CDF). Very roughly, what you know the fraction of the distribution $\Pr(y)$ that is below some $y$ for some set of values $\{y_i\}$. You can still estimate the log-normal parameters (or those of any other distribution) from this kind of information, but there are subtleties.

Two approaches come to mind. The first is a quick and dirty one that will not produce entirely accurate parameter estimates, but will get you close enough to get a sense of what the distribution looks like and, if you want, roughly what the Gini coefficient would be. The second is more complicated and more accurate for the kind of data you have.

Quick and dirty approximation

Here's the quick and dirty solution. The information you have is a "binned" version of the CDF, represented by a set of pairs $(q_i,y_i)$, where $q_i$ is the fraction of the distribution at or below the value $y_i$ (note: you said that the PPP is an average within the bin, which is technically distinct from the CDF, but for our calculation, that distinction doesn't make a difference).

Now, recall that the definition of the mean is $$\langle x \rangle = \sum_i x_i \Pr(x_i)\enspace ,$$ where $\Pr(x_i)$ is the probability of observing $x_i$. We don't have $\Pr(x)$, but we can approximate it using the binned CDF information, like this $$\hat{\mu} \approx \sum_{i=1}^k \Delta q_i \ln x_i$$ where $\Delta q_i=q_{i+1} -q_i$ is the size or width of the $i$th bin, out of $k$ bins. Similarly, for the standard deviation, the definition is $$\sigma = \sum_i (x_i-\langle x \rangle)^2 \Pr(x_i)\enspace,$$ which becomes $$\hat{\sigma} \approx \sum_{i=1}^{k} \Delta q_i (x_i-\hat{\mu})^2 \enspace .$$ To apply these to your data, you'll need to let $x_i=\ln y_i$ since you're working with the log-normal distribution, rather than the normal (or Gaussian) distribution. Coding up these estimators should be fairly easy.

In my numerical experiments with these estimators, I consistently get slight errors in the estimates relative to the underlying or "population" values I used to generate synthetic log-normal data. If you use these with your data, you should not treat the estimated values as being highly accurate. To get those, you'd need to apply a more mathematically sophisticated approach, which I'll sketch for you now.

Maximum likelihood approach

The more complicated and more accurate solution is to derive the maximum likelihood parameter estimate for the particular representation of the log-normal distribution you have, i.e., the binned CDF. The definition of the log-normal PDF is $$\Pr(x) = \frac{1}{x\sigma\sqrt{2\pi}}{\rm e}^{-\frac{(\ln x - \mu)^2}{2\sigma^2} } \enspace ,$$ and the CDF is $$\Pr(x<X) = F(x) = \frac{1}{2}\left(1+{\rm erf}\left( \frac{\ln x - \mu}{\sigma\sqrt{2}} \right) \right) \enspace ,$$ where $\textrm{erf}()$ is the error function, and where we let $F(x)$ be a short-hand representation for the CDF. (Normally, we would say $F(x\,|\,\mu,\sigma)$ to indicate that $F$ depends on your parameter choices, but I'm going to drop that notation henceforth; just remember that it's implied.) This is relevant because you want to assume your quantile data were drawn from a binned version of this distribution. If $F(x)$ is the CDF, i.e., the integral of $\Pr(x)$ from $-\infty$ to $x$, then let $F(x\,|\,a,b)$ be the integral of $\Pr(x)$ from $a$ to $b$. (Mathematically, $F(x\,|\,a,b)=F(b)-F(a)$.)

The log-likelihood of your observed quantile information is then $$\ln \mathcal{L} = \sum_{i=1}^k \ln F(x_i\,|\,q_i,q_{i+1})\enspace .$$ The more sophisticated approach would be to estimate $\mu$ and $\sigma$ by maximizing this function over these parameters. This would give you the maximum likelihood estimate for your log-normal model, given the observed information you have. For arbitrary choices of $\{q_i\}$, an analytic solution for the MLE is not possible, but for regularly spaced choices of the bin boundaries, it may be. Regardless, however, you may always numerically maximize the function (which many numerical software packages can do for you, if you whisper the right words to them).

What makes this approach more complicated is that you need to get the mathematics correct when you code up the numerical routine to do the estimation with the data. If the accuracy of your answers is really important, then this approach might be worth the extra effort.

-