7
$\begingroup$

I need to sample from the following mixture of two distributions:

$h_{\vec{\beta}}(r)=c(\vec{\beta})[(1-w_{m,\tau}(r))f_{\vec{\beta_{0}}}(r)+w_{m,\tau}(r)g_{\epsilon,\sigma}(r)]$

where $c(\vec{\beta})$ is a normalizing constant, $f_{\vec{\beta_{0}}}$ is the Gamma distribution, $g_{\epsilon,\sigma}$ is the Generalized Pareto distribution, and $w_{m,\tau}$ is the CDF of the Cauchy distribution: $1/2 +(1/\pi)*arctan[(r-m)/\tau]$, acting here as a weight/mixing function providing a smooth transition between the Gamma and Pareto (hence the adjective "dynamic" to qualify the mixture).

The weights allow having one unified mixture where the Gamma captures the bulk of the data, and as we get closer to the extremes, the Pareto takes over control. The approach originated in Figressi et al. (2002), and I have adopted the implementation presented in Vrac and Naveau (2007).

I have estimated the parameters $\vec{\beta}=(m,\tau,\gamma,\lambda,\epsilon,\xi)$, respectively the location and scale of the Cauchy CDF, shape and rate of the Gamma, and shape and scale of the GPD, following the procedure based on MLE described in the second paper referenced (R code can be found here), and now I want to sample from $h_{\vec{\beta}}$. Since the proportions of each component of the mixture are not fixed, I cannot just randomly sample from one distribution or the other (e.g., using using the rgamma() or rgpd() R functions) like is done here or here. So, how should I proceed?

If I could determine the CDF of the mixture, I could use inversion sampling ($U(0,1)$). EDIT: as noted by Xi'an in the second thread listed right below there is no closed form solution for the CDF.

EDIT: implementation of one of the two strategies in Xi'an's answer (rejection sampling) for my data here and Monte Carlo approach to estimating the quantile function of the mixture here.

$\endgroup$
  • 2
    $\begingroup$ Rejection sampling might work well here. If no one else comes along and does it first, I'll write up an answer later. $\endgroup$ – Dougal Jun 7 '15 at 18:30
  • $\begingroup$ Your density is well-defined (provided you can compute the normalising constant) and is a mixture in the usual sense but not of a Pareto and a Gamma. It writes as the sum of two terms, each can be renormalised as a density, and the normalising factors are the weights. But the two components of the mixture are non-standard. $\endgroup$ – Xi'an Jun 8 '15 at 11:36
8
$\begingroup$

As pointed out by Dougal, the shape of your target density$$h_\beta(r)\propto (1-w_{m,\tau}(r))f_{\beta_0}(r)+w_{m,\tau}(r) g_{\epsilon,\sigma}(r) $$is open to accept-reject simulation since $$(1-w_{m,\tau}(r))f_{\beta_0}(r)+w_{m,\tau}(r) g_{\epsilon,\sigma}(r) \le f_{\beta_0}(r)+g_{\epsilon,\sigma}(r)=2\left\{\frac{1}{2}f_{\beta_0}(r)+\frac{1}{2}g_{\epsilon,\sigma}(r)\right\} $$ Therefore simulating from the even mixture of Pareto $f_{\beta_0}$ and Gamma $g_{\epsilon,\sigma}$ and accepting with probability $$\dfrac{(1-w_{m,\tau}(r))f_{\beta_0}(r)+w_{m,\tau}(r) g_{\epsilon,\sigma}(r)}{f_{\beta_0}(r)+g_{\epsilon,\sigma}(r)}$$ would return you an exact output from your target density.

Note that the original paper by Frigessi et al. does include a way to simulate from the dynamic mixture on page 6: with probability $1/2$ simulate from $f_{\beta_0}$ and with probability $1/2$ from $g_{\epsilon,\sigma}$ [which is equivalent to simulating from the even mixture] and accept the outcome with probability $1-w_{m,\tau}(r)$ in the first case and $w_{m,\tau}(r)$ in the second case. It is unclear which one of those approaches has the highest average acceptance rate.

Here is a small experiment that shows the acceptance rates are comparable:

#Frigessi et al example
beta=2
lambda=gamma(1.5)
mu=tau=1
xi=.5
sigma=1
#the target is 
target=function(x) 
(1-pcauchy((x-mu)/tau))*dweibull(x,shape=beta,scale=1/lambda)+pcauchy((x-mu)/tau)*dgpd(x,xi=xi,beta=sigma)[1]

T=1e4
u=sample(c(0,1),T,rep=TRUE)
x=u*rweibull(T,shape=beta,scale=1/lambda)+(1-u)*rgpd(T,xi=xi,beta=sigma)
#AR 1
ace1=mean(runif(T)<(u*(1-pcauchy((x-mu)/tau))+(1-u)*pcauchy((x-mu)/tau)))
#AR 2
ace2=mean(runif(T)<target(x)/(dweibull(x,shape=beta,scale=1/lambda)+dgpd(x,xi=xi,beta=sigma)[1]))

with

> ace1
[1] 0.5173
> ace2
[1] 0.5473

An alternative is to use a Metropolis-Hastings algorithm. For instance, at each iteration of the Markov chain,

  1. pick the Pareto against the Gamma components with probabilities $1-w_{m,\tau}(x^{t-1})$ and $w_{m,\tau}(x^{t-1})$;
  2. Generate a value $y$ from the chosen component;
  3. Accept the value $y$ as $x^t=y$ with probability $$\dfrac{(1-w_{m,\tau}(y))f_{\beta_0}(y)+w_{m,\tau}(y) g_{\epsilon,\sigma}(y)}{(1-w_{m,\tau}(x^{t-1}))f_{\beta_0}(x^{t-1})+w_{m,\tau}(x^{t-1}) g_{\epsilon,\sigma}(x^{t-1})}$$ $$\times\dfrac{(1-w_{m,\tau}(y))f_{\beta_0}(x^{t-1})+w_{m,\tau}(y) g_{\epsilon,\sigma}(x^{t-1})}{(1-w_{m,\tau}(x^{t-1}))f_{\beta_0}(y)+w_{m,\tau}(x^{t-1}) g_{\epsilon,\sigma}(y)}$$ otherwise take $x^t=x^{t-1}$

The corresponding R code is straightforward

#MCMC style
propose=function(x,y){
#moving from x to y
target(y)*(pcauchy((y-mu)/tau,lowe=FALSE)*dweibull(x,shape=beta,scale=1/lambda)+pcauchy((y-mu)/tau)*dgpd(x,xi=xi,beta=sigma)[1:length(x)])/
(target(x)*(pcauchy((x-mu)/tau,lowe=FALSE)*dweibull(y,shape=beta,scale=1/lambda)+pcauchy((x-mu)/tau)*dgpd(y,xi=xi,beta=sigma)[1:length(x)]))}

x=seq(rgpd(1,xi=xi,beta=sigma),T)
for (t in 2:T){
  #proposal
  x[t]=rweibull(1,shape=beta,scale=1/lambda)
  if (runif(1)<pcauchy((x[t-1]-mu)/tau)) x[t]=rgpd(1,xi=xi,beta=sigma)
  #acceptance
  if (runif(1)>propose(x[t-1],x[t])) x[t]=x[t-1]}
ace3=length(unique(x))/T

and gives a higher acceptance rate

> ace3
[1] 0.877

While the fit is identical to the density estimate obtained by accept-reject: enter image description here

[Red curve for the accept-reject sample and blue curve for the MCMC sample, both based on 10⁴ original simulations]

$\endgroup$
  • $\begingroup$ Superb answer, many thanks. Sorry if I'm being slow, but regarding the rejection method, how do you determine the acceptance probability ($1/M$ using this notation). I know the probability you wrote is the area under the envelope (or instrumental) distribution $g(x)$ divided by the area under the arbitrary distribution $f(x)$, so intuitively it makes sense, but mathematically we just know that $f(x)<M*g(x)$ where $M>1$ and $P[u<f(x)/(M*g(x))]=1/M$ where $u$~$U(0,1)$ $\endgroup$ – Antoine Jun 10 '15 at 9:25
  • $\begingroup$ This is the beauty of accept-reject: if you can bound $f(x)<Mg(x)$, you do not need to separate $M$ from $g$ provided you can generate from $g$. This is exactly what happens in the second equation of my answer. I bounded the density of interest with 2 times a mixture with weights 1/2 and 1/2 that I can simulate. And I do not need the normalising constant of the lhs to do so. $\endgroup$ – Xi'an Jun 10 '15 at 9:48
  • 2
    $\begingroup$ This is definitely very sleek. I realized I was wrong in my first comment: the acceptance probability is not about areas, but simply values, or lengths. For each x sampled, if we draw a vertical line at x, the acceptance probability is the value of the arbitrary distribution over the value of the envelope distribution, which is the length of the segment $[0,f(x)]$ (= $f(x)$) divided by the length of the segment $[0,g(x)]$ (=$g(x)$). This is how you came up with the probability described in your third equation, right? $\endgroup$ – Antoine Jun 10 '15 at 10:21
  • $\begingroup$ Then, we simply uniformly sample along the vertical line and accept or reject with this probability... This is indeed very elegant. $\endgroup$ – Antoine Jun 10 '15 at 10:23
  • $\begingroup$ Yes, this is essentially the intuition behind accept-reject! $\endgroup$ – Xi'an Jun 10 '15 at 10:54
1
$\begingroup$

I suggest that you investigate Probabilistic Programming. Basically, you write programs that operate on probability distributions rather than (random) variables. This makes it easy to define functions that are dynamic mixtures of distributions. While full-featured Probabilistic Programming systems support Bayesian statistical inference, your application doesn't really need it. Primarily you just need the programming language and the Markov Chain Monte Carlo (MCMC) solvers to generate samples from the distribution defined by your program.

Good introductions and list of resources can be found in the following:

Good tutorials on Probabilistic Programming can be found at the Church web site, especially the second, "Simple Generative Models":

Finally, investigate STAN as an implementation platform. There are versions for R and Python. Again, you just need the generative model capabilities.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.