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!!

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

> qmixnorm(.95, c(0,2,5,10), c(1,1,1,1), c(.5,.25,.2,.05))
[1] 7.745521

• Unfortunately, this seems buggy at the moment, in terms of giving us a cryptic "collapsing to unique'x' values" error in KScorrect version 1.4.0. I have left a bug report. – Stephan Kolassa May 5 '20 at 10:00

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.

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