The answer to the long version with background:
This answer to the long version somewhat addresses another issue and,
since we seem to have difficulties formulating the model and the problem, I choose to rephrase it here, hopefully correctly.
For $1\le i\le I$, the goal is to simulate vectors $y^i=(y^i_1,\ldots,y^i_K)$ such that, conditional on a covariate $x^i$,
$$
y_k^i = \begin{cases} 0 &\text{ with probability }\operatorname{logit}^{-1}\left( \alpha_k x^i \right)\\
\log(\sigma_k z_k^i + \beta_k x^i) &\text{ with probability }1-\operatorname{logit}^{-1}\left( \alpha_k x^i \right)\end{cases}
$$
with $z^i=(z^i_1,\ldots,z^i_K)\sim\mathcal{N}_K(0,R)$. Hence, if one wants to simulate data from this model, one could proceed as follows:
For $1\le i\le I$,
- Generate $z^i=(z^i_1,\ldots,z^i_K)\sim\mathcal{N}_K(0,R)$
- Generate $u^i_1,\ldots,u^i_K \stackrel{\text{iid}}{\sim} \mathcal{U}(0,1)$
- Derive $y^i_k=\mathbb{I}\left\{u^i_k>\operatorname{logit}^{-1}\left( \alpha_k x^i \right)\right\}\, \log\{ \sigma_k z_k^i + \beta_k x^i\}$ for $1\le k\le K$
If one is interested in the generation from the posterior of $(\alpha,\beta,\mu,\sigma,R)$ given the $y^i_ k$, this is a harder problem, albeit feasible by Gibbs sampling or ABC.
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,
- pick the Pareto against the Gamma components with probabilities $1-w_{m,\tau}(x^{t-1})$ and $w_{m,\tau}(x^{t-1})$;
- Generate a value $y$ from the chosen component;
- 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:
[Red curve for the accept-reject sample and blue curve for the MCMC sample, both based on 10⁴ original simulations]
Best Answer
I'm assuming you are referecning to Inverse Transform Sampling method. Its very straight forward. Refer Wiki article and this site.
Pareto CDF is given by: $$ F(X) = 1 -(\frac{k}{x})^\gamma; x\ge k>0 \ and \ \gamma>0 $$
All you do is equate to uniform distribution and solve for $x$
$$ F(X) = U \\ U \sim Uniform(0,1) \\ 1 -(\frac{k}{x})^\gamma = u \\ x = k(1-u)^{-1/\gamma} $$
Now in R:
This will give you the random variate for Pareto distribution. I'm not sure about where you are getting
gammaroot?