Compute quantile function from a mixture of Normal distribution I have this mixture of normal distribution: 
$$X \sim \frac{1}{2}\mathcal{N}(\mu_{x_1}=10,\,\sigma_{x_1}^{2}=1)+\frac{1}{2}\mathcal{N}(\mu_{x_2}=13,\,\sigma_{x_2}^{2}=1)$$
How can i compute the quantile function of this distributionin R? 
I have to prove this function:
 qmist2n=function(y,mu1=10,sigma1=1,pi=0.5,mu2=13,sigma2=1){
  (1 - pi) * qnorm(y, mu1, sigma1) + pi * qnorm(y, mu2, sigma2)
}#qmist2n
I'm not sure if this code is correct.  Thanks for the help in advance!!
 A: Christoph Hanck's answer solves your problem and the link to Macro's answer explains it. 
But here is an additional answer to explain a bit more specifically about the error in your formula. 

With a mixture distribution you get indeed the situation that you can add upp the cumulative distribution function (see Macro's answer for the background).
For a mixture distribution $f_{Y}(x) = \sum a_i f_{X_i}(x)$ you have
$$F_{mixture}(x) = \sum a_i F_{X_i}(x)$$
but you do not have that the quantile function (the inverse $Q(x)=F^{-1}(x)$ of the cumulative distribution function) adds up similarly:
$$Q_{mixture}(x) \neq \sum a_i Q_{X_i}(x)$$
So your code is not correct because it is using a wrong formula.

Below is a graphic that may help to gain an intuitive idea about it:

You could see the graph of the quantile functions a sort of like a rotated image, flipping the x and y axis, of the graph of the cumulative distribution function. So you have that the two curves add up in the horizontal direction which is correctly done in the left graph (you add upp the curves for a given fixed $x$, I have added gray lines to highlight this, note that the lines are equal length on both sides, this is the case for equal weights 0.5 and 0.5). If you take a regular weighted mean of the quantile functions (combine the curves for a given quantile $q$ in the vertical direction), then you would get a sort of shift of the mean as you see in the right graph.

In R your could compute:
pmist2n = function(y,mu1=10,sigma1=1,pi=0.5,mu2=13,sigma2=1){
  (1 - pi) * pnorm(y, mu1, sigma1) + pi * pnorm(y, mu2, sigma2)
}#pmist2n

qmist2n = function(q,mu1=10,sigma1=1,pi=0.5,mu2=13,sigma2=1){
  # the min minmax below is computed to supply a range to the solver
  # the solution must be between the min and max
  # quantile of the mixed distributions
  minmax <- range(qnorm(q,mu1,sigma1),qnorm(q,mu2,sigma2)) 
  uniroot(function(x) pmist2n(x,mu1,sigma1,pi,mu2,sigma2)-q,
          interval = minmax,
          tol = 10^{-16})$root  
}#qmist2n

# example (if you look at the graph this must be below 10)
> qmist2n(0.2)
[1] 9.745184

This is a different approach from the KScorrect:qmixnorm function, which seems to approximate the quantile function by fitting a spline to points generated from the cdf (and somehow needs to generate 10000 random points in the process, times the number of normaldistributions in the mixture, but only uses the min and max of the distribution). This is not a very fast method, at least not on my machine (and using the version 1.2.0 of KScorrect).
benchmark("qmist2n" = {
            qmist2n(0.2)
          }, 
          "qmixnorm" = {
            qmixnorm(0.2,c(10,13),c(1,1),c(0.5,0.5))
          })
      test replications elapsed relative user.self sys.self user.child sys.child
1  qmist2n          100   0.041    1.000     0.040        0          0         0
2 qmixnorm          100  16.008  390.439    16.008        0          0         0

A: The question has already been discussed here.
Possibly after it's been dealt with in 2011, the following canned solution has become available in R:
library(KScorrect)
qmixnorm(.95, c(0,2,5,10), c(1,1,1,1), c(.5,.25,.2,.05))

Here, the first entry is the quantile, followed by the vector of means, standard deviations and weights (see https://rdrr.io/cran/KScorrect/man/dmixnorm.html).
The answer is as in the link:
> qmixnorm(.95, c(0,2,5,10), c(1,1,1,1), c(.5,.25,.2,.05))
[1] 7.745521

A: In general, you can compute the quantile function $Q(p) := \inf\{x:F(x)\geq p\}$ for any cdf $F(x)$ down to machine precision using the following approach:
# Computes the quantile function Q(p) for F(x) such that Q(p) = inf{x : F(x)>=p}
# Note that this is quantile(..., type=1) in R parlance
# It is assumed that exceedance probabilities are the first argument of pfun, other arguments are given in the ellipsis
quantile_bisect <- function(pfun, probs, ...) {
  
  opts <- list(...)
  
  probs[probs <= 0 | probs > 1] <- NA # Q is undefined in those cases
  if (all(is.na(probs))) {
    return(probs)
  }
  
  # sort probabilities for more efficient search, as we can reuse the lower bound L for larger p
  ord <- order(probs, na.last = T)
  probs <- probs[ord]
  
  # find initial range
  p <- probs[1]
  w <- 1 # search window size
  L <- 0
  while (do.call(pfun, c(L, opts)) >= p) {
    message(L)
    L <- L - w
    w <- 2 * w
  }
  R <- L
  
  res <- c()
  for (p in probs) {
    if (is.na(p)) {
      res <- c(res, NA)
    } else {
      while (do.call(pfun, c(R, opts)) < p) {
        L <- R
        R <- R + w
        w <- 2 * w # increase search window exponentially
      }
      # binary search (bisection) for p-quantile X, down to machine precision
      X = (L + R) / 2
      while (L < X & X < R) {
  
        if (do.call(pfun, c(X, opts)) >= p) {
          R = X
        } else {
          L = X
        }
        X = (L + R) / 2
      }
      res <- c(res, R)
    }
  }
  return(res[order(ord)])
}    


Fun <- function(x) {
  return((ecdf(1:10))(x))
}

probs <- seq(0, 1, length.out = 21)[2:21] # ecdf silently drops p=0, so we cannot compare the results
q <-  quantile(ecdf(1:10), probs, type = 1, na.rm = F) 
q_bisect <- quantile_bisect(Fun, probs)
all(q == q_bisect) # yes, we are comparing floats for strict equality here :-)

In terms of benchmarking, this might not be the most efficient solution. But as you can see in the last line, this approach also handles discontinuities and plateaus correctly, as opposed to many solutions that rely on defining a tolerance parameter. As floats are a finite set of rational numbers, not reals, the solutions are uniquely defined and computable, at least in theory (I'm sure there are a few numerical issues with this bisection approach as well, but you get the idea).
