# Get quantile function of dynamic mixture model

I have a dynamic mixture distribution fitted to my risk data (i.e., I have all parameters) of Weibull and Generalized Pareto, with a Cauchy CDF mixing function, that can be written as:


in R, this gives:

mixture = function(x){
((1 - pcauchy(x, location=535, scale=4.21e-04))*dweibull(x, shape=1.22, scale=62.31) +
pcauchy(x, location=535, scale=4.21e-04)*dgpd(x, xi=0.23, mu=0, beta=92.25))[1]
}


I want to know if the following is the correct way of writing the quantile function of my mixture (where the $q$ subscripts stand for the quantile functions):

\begin{align} &\mixture_{q}(y): y\in [0,1] \to \mixture_{q}(y) \in\mathbb{R^{+}} \\ &\mixture_{q}(y)=\big\{1-\Cauchy_{CDF}(\Weibull_{q}(y))\big\}\times \Weibull_{q}(y)+ \\ &\hspace{32mm} \Cauchy_{CDF}(GPD_{q}(y))\times GPD_{q}(y) \end{align}

in R, this translates to:

mixture.quantile = function(y){
((1 - pcauchy(qweibull(y, shape=1.22, scale=62.31), location=535, scale=4.21e-04)) *
qweibull(y, shape=1.22, scale=62.31) +
pcauchy(qgpd(y, xi=0.23, mu=0, beta=92.25), location=535, scale=4.21e-04) *
qgpd(y,xi=0.23,mu=0,beta=92.25))[1]
}


The only thing that changes is that because we are going from $[0,1]$ to $\mathbb{R^{+}}$, but $\Cauchy_{CDF}$ goes from $\mathbb{R^{+}}$ to $[0,1]$, the weighing has to be performed based on the values of the quantile functions and cannot be done based on the values of $x$ anymore (hence the $\Weibull_{q}(y)$ inside the $\Cauchy_{CDF}(.)$ for instance)... Is that correct??

With this quantile function I can get a Q-Q plot, or directly compare the quantile function of the mixture model with that of my data:

# load packages
libs = c("repmis","fExtremes","evmix","evir")
lapply(libs, library, character.only=TRUE)

data = source_data("https://www.dropbox.com/s/r7i0ctl1czy481d/test.csv?dl=0")[,1]

# get quantile function values
mixture.quant = numeric(length=length(data))
vector.quant  = ppoints(length(data), a=0)
for (i in 1:length(data)){
y                = vector.quant[i]
mixture.quant[i] = mixture.quantile(y)
}

# Q-Q plot
plot(mixture.quant, sort(data), xlab="theoretical quantile", ylab="sample quantile")
abline(sort(data), sort(data), col="red")

# empirical CDF with mixture quantile function overlaid
plot(mixture.quant, ppoints(length(data),a=0), col="red", type="l")
lines(quantile(data, ppoints(length(data),a=0)), ppoints(length(data),a=0), type="l")


PS: see these other related threads for background: 1, 2.

This is not correct: the cdf associated with the density $$(1-w_{\mu,\tau}(x))f_{\beta,\lambda}(x)+w_{\mu,\tau}(x)g_{\epsilon,\sigma}(x)$$is$$F(x)=\int_0^x \{(1-w_{\mu,\tau}(y))f_{\beta,\lambda}(y)+w_{\mu,\tau}(y)g_{\epsilon,\sigma}(y)\}\text{d}y$$hence is not the weighted sum of the cdfs of the Weibull and the GPD distributions. The same issue applies to the quantile function: it is not the weighted sum of the quantile functions of the Weibull and the GPD distributions. There is no closed form solution for either the pdf or the quantile functions, so the quantile has to be found numerically or by Monte Carlo (simulation).

• in your first sentence, I think you mean "the CDF associated with..." – Antoine Jun 18 '15 at 8:02
• so, the quantile function has to be solved numerically, that's the only way? – Antoine Jun 18 '15 at 8:05
• please see my "answer-comment" below: – Antoine Jun 18 '15 at 14:31

Note: this was too long to fit as a comment to Xi'an's answer, so I am posting it as an answer.

I have about 50,000 synthetic values generated from my mixture model using rejection sampling as described here:

# load packages
library(repmis)

sim=round(source_data("https://dl.dropboxusercontent.com/u/47464062/simulation.from.mixture.csv")[,1],2)


I can get the empirical CDF of these simulated values, along with its inverse, which is nothing else than the quantile function I'm looking for

ecdf.sim=ecdf(sim)
ecdf.sim(550)
quantile(sim,ecdf.sim(550))


Finally:

 # load historical data
data=source_data("https://www.dropbox.com/s/r7i0ctl1czy481d/test.csv?dl=0")[,1]

# get quantile function values
mixture.quant=numeric(length=length(data))
vector.quant=ppoints(length(data),a=0)
for (i in 1:length(data)){
y=vector.quant[i]
mixture.quant[i]=quantile(sim,y) # this is the only thing that changes
}


Using the same commands as in my question above, I can plot a Q-Q plot and the empirical CDF with the mixture CDF overlaid (not shown here for brevity). (Note: these plots are very close from the ones I obtained initially with the mathematically incorrect approach, but indicate a slightly better fit).

Xi'an, would this empirical way to tackle a parametric problem work? Can the resulting quantile function values be considered a good approximation of the true function?

• Correct: This is a Monte Carlo approach to evaluating the cdf and hence the quantile function of this distribution. You can compute the error in evaluating the cdf since this is a Binomial approximation, with variance p(1-p)/N, N being the number of simulations, hence plot an error band around the empirical cdf and re-hence invert this band into an error over the quantile. – Xi'an Jun 18 '15 at 17:01
• Now, your remark that the Q-Q plot is poor does not cover the same issue: you fit a dynamic mixture to your data and compare the empirical cdf based on the data with the (Monte Carlo version of the) theoretical cdf with estimated parameters. Since the data may (is most likely to) be from another distribution, there is no guarantee that the Q-Q plot shows a good fit. – Xi'an Jun 18 '15 at 17:03
• Since the true underlying distribution is unknown, would you recommend a fully nonparametric approach, with a Kernel Density Estimator for instance? (e.g., this thread of mine). – Antoine Jun 18 '15 at 20:32
• Using a non-parametric approach provides a fit to the data but lacks the explanatory features of a parametric model like the dynamic mixture. It all depends on the usage(s) you want to make of this estimate. – Xi'an Jun 18 '15 at 21:07
• you made a great point, this reminds me of this excellent paper by Leo Breiman. Ideally, being able to say: my data have been generated by this very model is very powerful. However, as you noted, and since I have very little knowledge about the underlying physical process at play, it is very likely that the data has not been generated by the formal model specified (especially with the quite bad fit that I'm getting). So, more humbly, it might be better not assuming anything and deal directly with the data only... – Antoine Jun 18 '15 at 21:14