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