Inverse CDF sampling for a mixed distribution

Question

The out-of-context short version

Let $y$ be a random variable with CDF $$ F(\cdot) \equiv \cases{\theta & y = 0 \\ \theta + (1-\theta) \times \text{CDF}_{\text{log-normal}}(\cdot; \mu, \sigma) & y > 0} $$

Let's say I wanted to simulate draws of $y$ using the inverse CDF method. Is that possible? This function doesn't exactly have an inverse. Then again there's Inverse transformation sampling for mixture distribution of two normal distributions which suggests that there is a known way to apply inverse transformation sampling here.

I'm aware of the two-step method, but I don't know how to apply it to my situation (see below).

---

The long version with background

I fitted the following model for a vector-valued response, $y^i = \left( y_1 , \dots , y_K \right)^i$, using MCMC (specifically, Stan):

$$ \theta_k^i \equiv \operatorname{logit}^{-1}\left( \alpha_k x^i \right), \quad \mu_k^i \equiv \beta_k x^i - \frac{ \sigma^2_k }{ 2 } \\ F(\cdot) \equiv \cases{\theta & y = 0 \\ \theta + (1-\theta) \times \text{CDF}_{\text{log-normal}}(\cdot; \mu, \sigma) & y > 0} \\ u_k \equiv F(y_k), \quad z_k \equiv\Phi^{-1}{\left( u_k \right)} \\ z \sim \mathcal{N}(\mathbf{0}, R) \times \prod_k f(y_k) \\ \left( \alpha, \beta, \sigma, R \right) \sim \text{priors} $$

where $i$ indexes $N$ observations, $R$ is a correlation matrix, and $x$ is a vector of predictors/regressors/features.

That is, my model is a regression model in which the conditional distribution of the response is assumed to be a Gaussian copula with zero-inflated log-normal marginals. I've posted about this model before; it turns out that Song, Li, and Yuan (2009, gated) have developed it and they call it a vector GLM, or VGLM. The following is their specification as close to verbatim as I could get it: $$ f(\mathbf{y}; \mathbf{\mu}, \mathbf{\varphi}, \Gamma) = c\{ G_1(y_1), \dots, G_m(y_m) | \Gamma \} \prod_{i=1}^m g(y_i; \mu_i, \varphi_i) \\ c(\mathbf{u} | \Gamma) = \left| \Gamma \right|^{-1/2}\exp\left( \frac{1}{2} \mathbf{q}^T \left( I_m - \Gamma^{-1} \right) \mathbf{q} \right) \\ \mathbf{q} = \left( q_1, \dots, q_m \right)^T, \quad q_i = \Phi^{-1}(u_i) $$ My $F_K$ corresponds to their $G_m$, my $z$ corresponds to their $\mathbf{q}$, and my $R$ corresponds to their $\Gamma$; the details are on page 62 (page 3 of the PDF file) but they're otherwise identical to what I wrote here.

The zero-inflated part roughly follows the specification of Liu and Chan (2010, ungated).

Now I would like to simulate data from the estimated parameters, but I'm a little confused as to how to go about it. First I thought I could just simulate $y$ directly (in R code):

for (i in 1:N) {
    for (k in 1:K) {
        Y_hat <- rbinom(1, 1, 1 - theta[i, k])
        if (Y_hat == 1)
            Y_hat <- rlnorm(1, mu[i, k], sigma[k])
    }
}

which doesn't use $R$ at all. I'd like to try to actually use the correlation matrix I estimated.

My next idea was to take draws of $z$ and then convert them back to $y$. This also seems to coincide with the answers in Generating samples from Copula in R and Bivariate sampling for distribution expressed in Sklar's copula theorem?. But what the heck is my $F^{-1}$ here? Inverse transformation sampling for mixture distribution of two normal distributions makes it sound like this is possible, but I have no idea how to do it.

Answer 1

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$,

1. Generate $z^i=(z^i_1,\ldots,z^i_K)\sim\mathcal{N}_K(0,R)$
2. Generate $u^i_1,\ldots,u^i_K \stackrel{\text{iid}}{\sim} \mathcal{U}(0,1)$
3. 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.

Answer 2

The answer to the out-of-context short version:

"Inverting" a cdf that is not invertible in the mathematical sense (like your mixed distribution) is feasible, as described in most Monte Carlo textbooks. (Like ours, see Lemma 2.4.) If you define the generalised inverse $$ F^{-}(u) = \inf\left\{ x\in\mathbb{R};\ F(x)\ge u \right\} $$ then $$ X \sim F \text{ is equivalent to } X=F^{-}(U)\text{ when } U\sim\mathcal{U}(0,1)\,. $$ This means that, when $F(y)$ has a jump of $\theta$ at $y=0$, $F^{-}(u)=0$ for $u\le\theta$. In other words, if you draw a uniform $\mathcal{U}(0,1)$ and it ends up smaller than $\theta$, your generation of $X$ is $x=0$. Else, when $u>\theta$, you end up generating from the continuous part, namely the log-normal in your case. This means using a second uniform random generation, $v$, independent from the first uniform draw and setting $y=\exp(\mu+\sigma\Phi^{-1}(v))$ to obtain a log-normal generation.

This is almost what your R code

Y_hat <- rbinom(1, 1, theta[i, k]) if (Y_hat == 1) Y_hat <- rlnorm(1, mu[i, k], sigma[k])

is doing. You generate a Bernoulli with probability $\theta_k^i$ and if it is equal to $1$, you turn it into a log-normal. Since it is equal to 1 with probability $\theta_k^i$ you should instead turn it into a log-normal simulation when it is equal to zero instead, ending up with the modified R code:

Y_hat <- rbinom(1, 1, theta[i, k])
    if (Y_hat == 0)
        Y_hat <- rlnorm(1, mu[i, k], sigma[k])