# Estimating a quantile from an unknown distribution

I have a variable $$x \sim N(0,1)$$, of dimension $$p$$, with each $$x_i$$ independent. I want to find the quantile, $$q$$, such that

$$\mathbb{P}\left(\frac{\sum_{i=1}^p (|x_i| - b_i) }{c} \geq q \right) = a,$$

where $$c>0$$ is a scalar, $$a \in (0,1)$$ is a specified value, and $$b$$ is a $$p$$-dimensional vector, such that $$b_1\geq \dots \geq b_p \geq 0$$.

As $$|x_i|$$ is a folded normal, I'm not sure how to derive an exact distribution, as I can't find any literature on the sum of folded normals. I believe we have a sum of independent folded-normals with unequal mean, but equal variance. I don't believe a closed-form exists, so then my question would be: how can I find $$q$$? I've tried to look into bootstrap, but without much success.

• Since the sum is taken over the vector components of vectors x and b and then divided by a scalar, the result is a scalar, too. What is the interpretation of the "greater equal" in comparison of the vector a? Nov 18, 2022 at 15:56
• @Dr_Be Sorry that was a mistake, a should be a scalar. Nov 18, 2022 at 15:58
• You are asking how to find the quantile function of a sum of iid Chi variates. There is no closed formula. There are good approximations for large $p.$
– whuber
Nov 18, 2022 at 23:15
• Cf. math.stackexchange.com/questions/2539306/…. From this I figure that you might not really have a chance to come up with an "exact" distribution, i.e. an explicit formula for the distribution (or its quantile function). Nov 19, 2022 at 11:22
• +1 for the quick response and improved question. You need numerical methods to find quantiles in this case.
– whuber
Nov 21, 2022 at 15:07

First let's simplify the question by defining $$q^* := cq + \sum_{i = 1}^pb_i$$ which then means we are asking for $$q^*$$ such that: $$P(\sum_{i = 1}^p|x_i| \geq q^*) = a$$

Now if you followed the link posted by @Dr_Be, than you mght have found this answer https://math.stackexchange.com/a/428790 , which sadly doesn't extend into higher dimnesions, as for $$p > 2$$, the set $$A := \{x \in \mathbb{R}^p : \sum_{i = 1}^p|x_i|< q\}$$ is a (Hyper-)Ocatahedron not a cube and so it can't be rotated as described. Now if i write out the integral exploiting a little bit of symmetry i get this:

$$P(\sum_{i = 1}^p|x_i| \geq q^*) = 1-2^p\int_{0}^{q^*}\int_{0}^{q^* - x_1}\cdots\int_{0}^{q^* - \sum_{i=1}^{p-1}x_i}\prod_{i = 1}^p\varphi(x_i) dx_p\cdots dx_1$$

You can gain a little bit traction there, but i don't see it being resolvable.

## Monte Carlo Estimation

So the exact answer will elude us and the integral looks very tough for most numerical solvers, but it is actually trivial to do with Monte Carlo simulations. Just simulate $$n$$ vectors $$x$$ and check for the empirical q. Here's an implementation in R:

n <- 10^5
p <- 20
a <- 0.8
z <- replicate(n, {
sum(abs(rnorm(p)))
})

q_star <- quantile(z, 1-a)


Quartiles are usually defines the other way around($$\leq$$), so i use 1-a. Now just recompute $$q = \frac{q^* - \sum_{i = 1}^pb_i }{c}$$

If you are interested in the standard error of this estimation i recommend the following chapter 4:A Tutorial on Quantile Estimation via Monte Carlo Hui Dong and Marvin K. Nakayama

If your $$p$$ is very large you can also use the CLT with each $$|x_i|$$ being a Chi-distribution with k = 1, https://en.wikipedia.org/wiki/Chi_distribution

## Implied question

The text below your equation seems to imply that what you actually want to look at is $$P(\sum_{i = 1}^p|x_i| \geq q) = a$$ with $$x_i \sim \mathcal{N}(b_i, c^2)$$. This is definitely hopeless with integrals as there is no longer any symmetry, but it barely complicates the Monte Carlo approach. Just use

b <- 1:p # put in your actual values
c <- 2
z <- replicate(n, {
sum(abs(rnorm(p, mean = b, sd = c)))
})


The CLT no longer applies as $$|x_i|$$ is no longer iid which can not be fixed by subtracting $$b_i$$s after the folding, as the folding has different effects depending on $$b_i$$.