8

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.

amoeba
  • 93,463
  • 28
  • 275
  • 317
user603
  • 21,225
  • 3
  • 71
  • 135
  • According to the notation and terminology, we are supposed to understand "$L_z(t)$" as mapping any $F$ into the *number* $L_z(t)[F]=P_F(|z-Z|\le t)$ (that's what "functionals" do, after all). But then what could "$L_z^{-1}(q)$" possibly mean? It could be either a *distribution* (literally, the inverse of $L_z(t)$) or, more likely, a *number* (if we fix $F$ and view $L_z$ as a *function*--not functional!), but I don't see any way to interpret it as a *random variable* with the distribution $F$. That makes $M_q(t)$ fairly mysterious. – whuber Jan 16 '15 at 22:13
  • A tiny bit, but there must still be some typos. Perhaps in the definition of $L(t)$ you meant to use "$Z$" (a random variable with $F$ for its distribution) instead of "$z$"? After all, one cannot assign a probability to "$H(z)\le t$", since both $H(z)$ and $t$ are *numbers.* You are running a big risk here that a single tiny typo in any equation will turn your question into something completely unintended and irrelevant; even worse, you might get correct answers that *you* misinterpret! It would help (a lot) to include English explanations of what each formula is intended to represent. – whuber Jan 16 '15 at 22:27

3 Answers3

1

$\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.

Ben Kuhn
  • 5,373
  • 1
  • 16
  • 27
  • 1
    I don't think this addresses the problem, because we are not dealing with realizations of $F$ here but with $F$ itself (If you will, I m interested in the values of these objects when the sample size goes to $\infty$). But maybe I misunderstand your answer? – user603 Jan 17 '15 at 10:17
  • 1
    Well, I just realized I made a mistake with the derivation anyway--the optimizations are actually nested. I think there are still optimization techniques you can use, but I don't know if they're better than just taking the second median on a large sample as you've been doing. – Ben Kuhn Jan 17 '15 at 18:50
  • good so we now wonder the same thing;) – user603 Jan 17 '15 at 18:51
1

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.

Community
  • 1
Antoine
  • 5,740
  • 7
  • 29
  • 53
  • Why do you think you get better precision by bootstrapping? – kjetil b halvorsen Oct 08 '15 at 16:19
  • the fate of every finite sample is that it does not contain the full spectrum of observations that may occur. E.g., the max stream flow of a river observed over 100 years is obviously not the absolute max that can happen. So your estimates of the 500-year flood (0.998 quantile) or 1000-year flood (0.999 quantile) based on your limited sample will be biased (risk will be **underestimated**). On the opposite, if you generate hundreds of thousands of new observations by simulating (via smoothed bootstrap or any other technique), your estimates will be more accurate – Antoine Oct 09 '15 at 08:03
  • That is an misunderstanding! The values in the bootstrapped sample do all stem from that same limited, finite sample, and doesnt contain more information than the sample itself. Bootstrapping (like other analytical techniques) can only help us understand better what information is in the sample, it cannot increase that information. – kjetil b halvorsen Oct 09 '15 at 08:15
  • No! I am talking about the **smoothed** bootstrap. It generates never-seen observations, that exceed the range of the original sample. Please click the link in my answer above. – Antoine Oct 09 '15 at 08:16
  • @ Antoine: OK, but that does not change anything in my comment. Those smoothed "obs" are generated from your model of the phenomenon, not the phemonenon itself, so is not "data". – kjetil b halvorsen Oct 09 '15 at 08:18
  • `Those smoothed "obs" are generated from your model of the phenomenon` is wrong. The simulated values are generated from the **sample**, not some model. This is a fully data-driven, non-parametric approach. But in any case, your critic is true for any Monte Carlo technique. Of course, when you simulate you're not generating "true" observations, but it's better than nothing. If you had plenty obs, you wouldn't need to simulate in the first place. – Antoine Oct 09 '15 at 08:24
  • Well, your smoothing of the sample is based on some model, isn't it? – kjetil b halvorsen Oct 09 '15 at 08:26
  • No! Really, please take the time to go through the links I provided. I'm not going to repeat everything in comments. – Antoine Oct 09 '15 at 08:29
  • Why not leverage quantile regression? It is nonparametric, distribution-free, robust to extreme values and, using Rob Hyndman's method described in this paper, you can fit all of the quantiles in a grid-like fashion. Once you have the full predicted distribution, you can assign any underlying probabilistic distribution that makes sense to you. robjhyndman.com/papers/smart-meter-quantiles.pdf – Mike Hunter Oct 09 '15 at 09:57
  • Actually, I have managed the solve this one. I will post a full answer in a bit. – user603 Oct 09 '15 at 13:32
1

So, I think that the best way to obtain

$$\text{med}_{Z\sim F} H(Z)$$

is to:

  1. 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))$
  2. Compute the weighted median of $\{H(z_i)\}_{i=1}^n$ with weights $F_Z^\prime(z_i)$.
user603
  • 21,225
  • 3
  • 71
  • 135