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.

  • $\begingroup$ @Xi'an it's a Gaussian copula, for estimating the dependence among the $y$ components. $\endgroup$ Dec 16 '14 at 17:44
  • 1
    $\begingroup$ The thread you reference about sampling from mixtures of Normals applies directly to your problem with no essential modification: instead of using the inverse CDFs of Normals, use the inverse CDFs of your two components. The inverse CDF of the atom at $y=0$ is the constant function, always equal to $0$. $\endgroup$
    – whuber
    Dec 16 '14 at 18:24
  • $\begingroup$ @whuber I'm just confused about how to use the inverse CDFs of the two components: what do I draw, what do I draw it from, and then what do I plug each thing into? $\endgroup$ Dec 16 '14 at 19:07
  • 1
    $\begingroup$ @Xi'an nicely explains that in his answer to the normal-mixture question: you use a uniform variate to select the mixture component and then you draw a value from that component (in any way you like). In your case it's exceptionally easy to draw a value from the first component: it's always $0$! To draw a value from the second component use any lognormal random number generator you like. In each case you wind up with a number: there's no "plugging in" to accomplish; the whole objective of random number generation is to obtain that number. $\endgroup$
    – whuber
    Dec 16 '14 at 19:11
  • $\begingroup$ @whuber the new answer cleared it up for me. Thank you both. $\endgroup$ Dec 16 '14 at 19:14

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.

  • 1
    $\begingroup$ I knew I was missing something. "Everything is obvious in hindsight." My intent: I'm interested in the value of $F(y^i | x^i)$, so yes, I'm interested in drawing from the joint posterior of the parameters. I want the simulated $y$'s to see if the model fits. $\endgroup$ Dec 16 '14 at 19:12
  • 1
    $\begingroup$ How is the second problem much harder? I've already estimated the model and I have posterior draws. We can continue in chat if you want, to keep from cluttering up the comments here. $\endgroup$ Dec 16 '14 at 19:20
  • 1
    $\begingroup$ Oh, in general, yeah. Fortunately I have Stan and the No-U-Turn Sampler doing the hard work for me there. $\endgroup$ Dec 16 '14 at 20:23

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])
  • $\begingroup$ So all together, my simulation procedure would be: 1) draw $z$, 2) compute $u_k = \Phi(z_k)$, then 3) compute $y_k = 0$ if $u_k \leq \theta$ and $y_k = F^{-1}_\text{log-normal}(u_k)$ otherwise. Correct? $\endgroup$ Dec 16 '14 at 17:31
  • $\begingroup$ No, incorrect. You draw a first uniform to decide between $0$ and log-normal, then a second uniform in case you have decided for a log-normal. See the edited version of my answer. $\endgroup$
    – Xi'an
    Dec 16 '14 at 17:33
  • $\begingroup$ But that ignores the $z$ component; hence my question. I made a clarifying edit, and also addressed the mistake in my pseudocode. $\endgroup$ Dec 16 '14 at 17:44
  • $\begingroup$ My answer is for the short version and for the R code you provided. I hope it helps for the long-version, but your formula for the joint model is still incorrect. You should define the model on the $y$'s without using the uniforms... $\endgroup$
    – Xi'an
    Dec 16 '14 at 17:59
  • $\begingroup$ How is that model incorrect? I just plugged my $F_1, \dots, F_K$ into the formula provided by the paper I cited (corresponding to $G_1, \dots, G_m$ in their notation). Is that invalid? $\endgroup$ Dec 16 '14 at 18:51

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.