Fit a parametrized distribution on a set of quantiles

Question

first of all, I'm not a statistician nor a data scientist, but a software developer. Thus, although I do have some (old) knowledge in statistics and probabilities, my vocabulary may not be very accurate.

I'm working on an app in which we must draw incomes from a parametrized distribution such as lognormal, gamma, Weibull, etc. (the choice is to be made yet). We need to fit that distribution before we raw values. The available input data are quantiles, most often (but not always) the 1st to 9th deciles.

I'm looking for a mathematical way to perform the fit that wil later be implemented in Python/SciPy. I know there are least-square methods to do so but I'd prefer to use maximum-likelihood methods. What I was thinking of is the following. Let us call $q_i$ the available (observed) quantiles, each for a "rank" $x_i$. If the sample is $N$-sized, and with $F$ being the CDF of the distribution, then the probability to draw $Nx_i$ values less or equal to $q_i$ is: $$P(n=Nx_i)={{N}\choose{Nx_i}}F(q_i)^{Nx_i}(1-F(q_i))^{N(1-x_i)}$$ From there I was thinking to maximize the following "likelihood" estimator: $$\mathcal{L(a_1, a_2, ...)}=\sum_i x_i\ln{F(q_i; a_1, b_2, ...)}+(1-x_i)\ln{(1-F(q_i; a_1, a_2, ...))}$$ where $a_1$, $a_2$, $...$ are the distribution parameters.

Now, admittedly, I'm not quite sure the approach is correct. What do you think of it? In particular, should I consider not the likelihood for $Nx_i$ values to belong to $(-\infty, q_i]$ but rather the likelihood for $N(x_{i+1}-x_i)$ values to fall in the $(q_i, q_{i+1}]$ interval?

Also, feel free to correct my terms and considerations. E.g., I'd be glad to learn how to call the $x_i$'s rigorously.

EDIT: here comes another estimator to maximize given $I$ observed quantiles $q_0$, $q_1$... $q_{I-1}$. This one is based on the likelihood to find $Nx_0$ items up to $q_0$, then $N(x_1-x_0)$ of the remaining items in the $(q_0, q_1]$ interval, and so on: $$\mathcal{L}_2=x_0\ln{F(q_0)}+(1-x_{I-1})\ln(1-F(q_{I-1}))+\sum_{i=0}^{i=I-2}(x_{i+1}-x_i)\ln(F(q_{i+1})-F(q_i))$$ Would $\mathcal{L}_2$ be more or less relevant than $\mathcal{L}$? Explanations will obviously be welcome!

Answer 1

Thanks to whuber for their suggestion! Reading the suggested post, I thought of the following estimator to look consistent and it appears to work pretty well. It is based on the number of sample values expected to be found in each quantile interval.

The problem to solve here is the following. You are given a set of $I$ quantiles $q_1$, $q_2$... $q_I$, all distinct, sorted by increasing values and corresponding to known cumulative probabilities $x_1$, $x_2$... $x_I$, all in the $(0, 1)$ interval. For example, if $x_5=0.5$, then $q_5$ is the median. You suspect that the quantiles can be fitted with some continuous distribution with parameters $a_1$, $a_2$, etc. What you need is to find the parameter values that will maximize the likelihood for the distribution to account for the quantiles provided.

Let us call $P$ the distribution's probability density function (PDF), and $F$ the cumulative density function (CDF). If the quantiles are extracted from an $N$-sized population, then $Nx_1$ values were found to be less than or equal to $q_1$ (well, in practice, $q_1$ was probably interpolated between two sample values). Likewise, $N(x_2-x_1)$ values were found in the $(q_1, q_2]$ interval, and so on. On the other hand, from $F$, you can derive the theoretical probability for a given sample item to belong to a given quantile interval. For example, that probability will be $F(q_2)-F(q_1)$ for interval $(q_1, q_2]$. From there, you can establish a multinomial probability for each scenario on how the quantile intervals will be populated.

For the sake of computation, let us add two cumulative probabilities to our list: $x_0=0$ and $x_{I+1}=1$. Whatever the distribution's parameters, we must force $F(x_0)=0$ and $F(x_1)=1$. Working with logarithms of probability/likelihood, we can come to the likelihood estimator: $$\mathcal{L}(a_1, a_2...)=\sum_{i=0}^{I}(x_{i+1}-x_i)\ln{\frac{F(q_{i+1}; a_1, a_2...)-F(q_i; a_1, a_2...)}{x_{i+1}-x_i}}$$ Note the division of $F(q_{i+1})-F(q_i)$ by $x_{i+1}-x_i$. The motive is that in the case of a perfect fit, $F(q_i)=x_i$ so $\mathcal{L}$ is exactly zero. With real-world data, expect $\mathcal{L}<0$. If $P$ is defined in the whole quantile range, then $F$ can be derivated, and so can $\mathcal{L}$. Hence, so as to maximize $\mathcal{L}$, we'll want to find parameter values so that the derivatives of $\mathcal{L}$ with respect to the parameters is zero. The derivative with respect to $a_j$ is: $$\frac{\partial \mathcal{L}}{\partial a_j}(a_1, a_2...)=\sum_{i=0}^{I}\frac{x_{i+1}-x_i}{F(q_{i+1}; a_1, a_2...)-F(q_i; a_1, a_2...)}\left(\frac{\partial F(q_{i+1}; a_1, a_2...)}{\partial a_j}-\frac{\partial F(q_{i}; a_1, a_2...)}{\partial a_j}\right)$$ with the partial derivatives forced to zero at $i=0$ and $i=I+1$.

Finally, let us have a look to the case where the parameters are location $\mu$ and scale $s$. Calling $P_0$ and $F_0$ the PDF and CDF at $\mu=0$ and $s=1$, we have: $$P(q; \mu, s)=\frac{1}{s}P_0\left(\frac{q-\mu}{s}\right)$$ $$F(q; \mu, s)=F_0\left(\frac{q-\mu}{s}\right)$$ Hence, the partial derivatives of $F$ with respect to $\mu$ and $s$ are: $$\frac{\partial F}{\partial \mu}(\mu, s)=-\frac{1}{s}P_0\left(\frac{q-\mu}{s}\right)=-P(q; \mu, s)$$ $$\frac{\partial F}{\partial s}(\mu, s)=-\frac{1}{s}\frac{q-\mu}{s}P_0\left(\frac{q-\mu}{s}\right)=-\frac{q-\mu}{s}P(q; \mu, s)$$ Hope this is useful to somebody else...