# How to estimate the third quartile of binned data?

Is there any technical trick to determine the third quartile if it belongs to an open interval which contains more that one fourth of the population (so I can't close the interval and use the standard formula)?

### Edit

In case I misunderstood something I will provide more or less full context. I have data arranged in a table with two columns and, say, 6 rows. With each column corresponds an interval (in the first column) and a quantity of population which "belongs" to that interval. The last interval is open and includes more than 25% of the population. All intervals (with exception of the last) have the same range.

Sample data (transposed for presentation):

Column 1: (6;8),(8;10),(10;12),(12;14),(14;16),(16;∞)
Column 2:    51,    65,     68,     82,     78,   182


The first column is to be interpreted as an income level range. The second is to be interpreted as the number of employees whose income belongs to the interval.

The standard formula I'm thinking about is $$\mathbb{Q}_{3}=x_{Q_{3}}+ \frac{\frac{3N}{4}- \sum_{i=1}^{k-1}n_{i}}{n_{Q_{3}}}r_{Q_{3}}$$.

• A common assumption when trying to estimate quantiles with binned data is to assume uniformity within bins. But when you know something about the way the data is likely to be distributed (as with incomes, which are right skew) assumptions that reflect that knowledge will tend to be better. Another alternative would be to assume that it's smooth, and then smooth the data (whether by KDE or some fitted distribution), redistribute points within bins according to the model [& possibly re-estimate (in somewhat EM-like fashion) the fit, &redistribute in bins again] then estimate quantiles from that. – Glen_b Apr 27 '14 at 3:48

You need to fit these binned data with some distributional model, for that is the only way to extrapolate into the upper quartile.

### A model

By definition, such a model is given by a cadlag function $F$ rising from $0$ to $1$. The probability it assigns to any interval $(a,b]$ is $F(b)-F(a)$. To make the fit, you need to posit a family of possible functions indexed by a (vector) parameter $\theta$, $\{F_\theta\}$. Assuming that the sample summarizes a collection of people chosen randomly and independently from a population described by some specific (but unknown) $F_\theta$, the probability of the sample (or likelihood, $L$) is the product of the individual probabilities. In the example, it would equal

$$L(\theta) = (F_\theta(8) - F_\theta(6))^{51} (F_\theta(10) - F_\theta(8))^{65} \cdots (F_\theta(\infty) - F_\theta(16))^{182}$$

because $51$ of the people have associated probabilities $F_\theta(8) - F_\theta(6)$, $65$ have probabilities $F_\theta(10) - F_\theta(8)$, and so on.

### Fitting the model to the data

The Maximum Likelihood estimate of $\theta$ is a value which maximizes $L$ (or, equivalently, the logarithm of $L$).

Income distributions are often modeled by lognormal distributions (see, for example, http://gdrs.sourceforge.net/docs/PoleStar_TechNote_4.pdf). Writing $\theta = (\mu,\sigma)$, the family of lognormal distributions is

$$F_{(\mu, \sigma)}(x) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{(\log(x)-\mu)/\sigma} \exp(-t^2/2) dt.$$

For this family (and many others) it is straightforward to optimize $L$ numerically. For instance, in R we would write a function to compute $\log(L(\theta))$ and then optimize it, because the maximum of $\log(L)$ coincides with the maximum of $L$ itself and (usually) $\log(L)$ is simpler to calculate and numerically more stable to work with:

logL <- function(thresh, pop, mu, sigma) {
l <- function(x1, x2) ifelse(is.na(x2), 1, pnorm(log(x2), mean=mu, sd=sigma))
- pnorm(log(x1), mean=mu, sd=sigma)
logl <- function(n, x1, x2)  n * log(l(x1, x2))
sum(mapply(logl, pop, thresh, c(thresh[-1], NA)))
}

thresh <- c(6,8,10,12,14,16)
pop <- c(51,65,68,82,78,182)
fit <- optim(c(0,1), function(theta) -logL(thresh, pop, theta[1], theta[2]))