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.

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