How to obtain the quantile function when an analytical form of the distribution is not known The problem comes from page 377-379 of this [0] paper. 
Given a continuous distribution $F$ and a fixed $z\in\mathbb{R}$,
 consider:
$$L_z(t)=P_F(|z-Z|\leq t)$$
and 
$$H(z)=L^{-1}_z(0.5)=\underset{Z\sim F}{\mbox{med}}|z-Z|$$
where $L^{-1}_z(u)=\inf\{t:L_z(t)>u\}$ is the right continuous inverse.
So for a fixed $z$, this is the median distance of all 
$Z\sim F$ to $z$. Next, consider the function:
$$L(t)=P_F(H(Z)\leq t)$$ 
Now, I don't have an analytical expression for $H(z)$ (in fact I'm pretty sure an analytical expression for it is not possible) but given a CDF $F$ I can easily uses a root finding algorithm to obtain $H(z)$ for any given $z$.
In this application, the interest is on:
$$L^{-1}(0.5)=\underset{Z\sim F}{\mbox{med}}H(Z)$$
This is the median value of the $H(Z)$, again, for $Z\sim F$.
Right now to get $L^{-1}(0.5)$, I compute (as explained above, using a root finding algorithm) values of $H(z)$ corresponding to many values of $z$ on a grid and take the weighted median of these values of $H(z)$ (with weights $f(z)$) as my estimate of $L^{-1}(0.5)$. 
My questions are: 


*

*Is there a more accurate approach to get $L^{-1}(0.5)$ (the  authors of the paper do not say how $L^{-1}(0.5)$ is computed) and

*How should the grid of values of $z$ be chosen?
[0] Ola Hössjer, Peter J. Rousseeuw and Christophe Croux.  Asymptotics of an estimator of a robust spread functional. Statistica Sinica 6(1996), 375-388.
 A: $\DeclareMathOperator*{\med}{med}$The median is the point that minimizes the expected $L^1$ distance:
$$\med_Z f(Z) = \arg\min_m E_z|f(Z) - m|$$
Hence we can simplify your expression:
$$\begin{equation}\med_{z_1 \sim F} \med_{z_2 \sim F} |z_1 - z_2| \\
= \arg\min_{m_1}E_{z_1 \sim F}\left| m_1 - \arg\min_{m_2} E_{z_2 \sim F}\left| m_2 - \left|z_1 - z_2\right|\right|\right|
\end{equation}$$
I think this is a bilevel optimization problem, which I don't know too much about but perhaps there are standard techniques you can apply. Then again, it might not be any faster than just calculating the sample median of medians for larger samples until convergence.
A: A straightforward data-driven approach to estimating the quantile function consists in:


*

*bootstrapping your observations to generate many more values than are
in your original sample (especially, values beyond the range of the initial limited sample). A good strategy is to use a smoothed bootstrap simulation
scheme to avoid the main limitations of the basic nonparametric
bootstrap. This is equivalent to simulating from a Kernel Density Estimate.

*from this, you can get the empirical Cumulative Distribution Function
(CDF) of the simulated values (ecdf function in R). The inverse of the CDF is nothing else than the quantile function (quantile
function in R). See here to get the values and plot your quantile
function. You can even get confidence bands.


A pre-requisite though is that you sample features enough observations to at least get a good idea of the shape of your underlying PDF.
A: So, I think that the best way to obtain 
$$\text{med}_{Z\sim F} H(Z)$$ 
is to:


*

*compute the entries of the $n$ vector $\{H(z_i)\}_{i=1}^n$ of values of $H(z_i)$ corresponding to a grid of $n$ values of $\{z_i\}_{i=1}^n$ placed uniformly on $(F_Z^{-1}(\epsilon),F_Z^{-1}(1-\epsilon))$

*Compute the weighted median of $\{H(z_i)\}_{i=1}^n$ with weights 
$F_Z^\prime(z_i)$.

