# Gaussian copula for discrete marginals

**URL:** <https://discourse.mc-stan.org/t/gaussian-copula-for-discrete-marginals/691>\
**Category:** Modeling\
**Created:** [May 29, 2017, 1:44pm UTC](https://discourse.mc-stan.org/t/gaussian-copula-for-discrete-marginals/691 "2017-05-29T13:44:41Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![Ben\_Lambert](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ben_lambert/32/4407_2.png) [@Ben\_Lambert](https://discourse.mc-stan.org/u/Ben_Lambert)\
**Post date:** [May 29, 2017, 1:44pm UTC](https://discourse.mc-stan.org/t/gaussian-copula-for-discrete-marginals/691/1 "2017-05-29T13:44:41Z")

</div>

Hi,

I would like to estimate a model that uses a Gaussian copula to piece together discrete marginals for count data (e.g. poisson or negative binomial distributed). To do this requires the multivariate Gaussian, which I don’t think currently exists in Stan. (Please correct me if I am wrong, though.)

It also requires the inverse normal cdf, which I also don’t think exists in Stan? Again, shoot me down if I am wrong.

I attach the paper that I am working off. In particular I would like to use eqn. (2.5) to work out the probability of a discrete data point.

If anyone has any bright ideas how to do this I am all ears! (I don’t want to resort to using multivariate poissons etc. – my problem has too many dimensions for this to be feasible.) Otherwise I will go and code up the Gibbs sampler in section 4 of the paper.

Best,

Ben

[Gaussian copula better.pdf](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/1X/d18c32d00e48705997e789ab7341480a26a1694a.pdf) (272.2 KB)

---

<div class="post-metadata">

**Author:** ![bgoodri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bgoodri/32/4451_2.png) [@bgoodri](https://discourse.mc-stan.org/u/bgoodri)\
**Post date:** [May 29, 2017, 2:51pm UTC](https://discourse.mc-stan.org/t/gaussian-copula-for-discrete-marginals/691/2 "2017-05-29T14:51:32Z")

</div>

> [@Ben\_Lambert](#):
>
> To do this requires the multivariate Gaussian, which I don’t think currently exists in Stan. (Please correct me if I am wrong, though.)

We don’t have the Gaussian copula density if that is what you are asking, but you can use `multi_normal_prec` with precision matrix equal to Gamma inverse minus the identity matrix and then correct the determinant afterwards. But you would probably be better off writing your own function that takes a Cholesky factor of a correlation matrix.

> [@Ben\_Lambert](#):
>
> It also requires the inverse normal cdf, which I also don’t think exists in Stan? Again, shoot me down if I am wrong.

This is totally wrong; it is the `inv_Phi` function.

---

<div class="post-metadata">

**Author:** ![Ben\_Lambert](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ben_lambert/32/4407_2.png) [@Ben\_Lambert](https://discourse.mc-stan.org/u/Ben_Lambert)\
**Post date:** [May 29, 2017, 2:56pm UTC](https://discourse.mc-stan.org/t/gaussian-copula-for-discrete-marginals/691/3 "2017-05-29T14:56:08Z")

</div>

Hi Ben,

Thank you, and my apologies for the misunderstanding here – I meant to say multivariate gaussian CDF, rather than just multivariate gaussian, (which I don’t think exists in Stan?).

Also apologies for not knowing where to look in the pages of the Stan manual for the univariate inverse normal cdf. I could not find it next to normal\_cdf and so assumed it wasn’t there.

Best,

Ben

---

<div class="post-metadata">

**Author:** ![bgoodri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bgoodri/32/4451_2.png) [@bgoodri](https://discourse.mc-stan.org/u/bgoodri)\
**Post date:** [May 29, 2017, 3:20pm UTC](https://discourse.mc-stan.org/t/gaussian-copula-for-discrete-marginals/691/4 "2017-05-29T15:20:49Z")

</div>

The multivariate Gaussian CDF does not exist anywhere really. But you don’t need it to do the stuff in section 4 of that paper.

---

<div class="post-metadata">

**Author:** ![Ben\_Lambert](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ben_lambert/32/4407_2.png) [@Ben\_Lambert](https://discourse.mc-stan.org/u/Ben_Lambert)\
**Post date:** [May 29, 2017, 3:53pm UTC](https://discourse.mc-stan.org/t/gaussian-copula-for-discrete-marginals/691/5 "2017-05-29T15:53:00Z")

</div>

Ok, thank you. Sorry am still being a bit slow though (first time working with Copulas) – to code this up in Stan would it be ok to work with with augmented likelihood (eqn. at the top of section 4.1), setting priors on Gamma and Theta?

Best, Ben

---

<div class="post-metadata">

**Author:** ![bgoodri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bgoodri/32/4451_2.png) [@bgoodri](https://discourse.mc-stan.org/u/bgoodri)\
**Post date:** [May 29, 2017, 4:27pm UTC](https://discourse.mc-stan.org/t/gaussian-copula-for-discrete-marginals/691/6 "2017-05-29T16:27:46Z")

</div>

You don’t want to do that indicator function stuff with NUTS. Rather, you want to declare a bunch of nuisance parameters that are between 0 and 1, and then adjust them to the relevant bounds, and then make the Jacobian adjustment. It is a bit like what I outlined here:

[https://groups.google.com/d/msg/stan-users/GuWUJogum1o/LvxjlUBnBwAJ](https://groups.google.com/d/msg/stan-users/GuWUJogum1o/LvxjlUBnBwAJ)

---

<div class="post-metadata">

**Author:** ![Ben\_Lambert](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ben_lambert/32/4407_2.png) [@Ben\_Lambert](https://discourse.mc-stan.org/u/Ben_Lambert)\
**Post date:** [May 30, 2017, 12:32am UTC](https://discourse.mc-stan.org/t/gaussian-copula-for-discrete-marginals/691/7 "2017-05-30T00:32:47Z")

</div>

Hi Ben,

Thank you very much for your response, and I must say I did enjoy your approach to sampling from a truncated multivariate normal in the other question.

I do see the similarities between this problem and the copula one, however I am still not quite joining all the dots. Are you suggesting to do something like the following,

- Create a load of u\_ij which are uniformly distributed between the following bounds: a\_ij \<= u\_ij \< b\_ij, where a\_ij = F\_j(y\_ij-,theta\_ij) and b\_ij = F(y\_ij,theta\_ij). Here F(y\_ij-) = F(y\_ij-1) since y\_ij is discrete. (Here F(.) is the CDF for a discrete distribution, say the poisson.)
- Transform these to x\_ij = Phi^{-1}(u\_ij).
- Calculate the likelihood for a N( **x** ,| **0** , **Gamma** ), and use this to increment the log probability.
- Correct for the non-linear transformation of u - \> x via a Jacobian.

I would be really grateful if you have the chance to look over the above, and tell me your thoughts. I do suspect that these types of model would be really helpful for people in Stan as they look for alternatives to multivariate discrete distributions.

Best,

Ben

---

<div class="post-metadata">

**Author:** ![Ben\_Lambert](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ben_lambert/32/4407_2.png) [@Ben\_Lambert](https://discourse.mc-stan.org/u/Ben_Lambert)\
**Post date:** [May 30, 2017, 6:25pm UTC](https://discourse.mc-stan.org/t/gaussian-copula-for-discrete-marginals/691/8 "2017-05-30T18:25:40Z")

</div>

Just following up on my own point here: would the likelihood be an unbounded N( **x** , **Gamma** ) distribution, or a truncated one? Or am I still way off here?

Best,

Ben

---

<div class="post-metadata">

**Author:** ![Ben\_Lambert](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ben_lambert/32/4407_2.png) [@Ben\_Lambert](https://discourse.mc-stan.org/u/Ben_Lambert)\
**Post date:** [May 31, 2017, 2:17am UTC](https://discourse.mc-stan.org/t/gaussian-copula-for-discrete-marginals/691/9 "2017-05-31T02:17:53Z")

</div>

So the following R code generates data from a bivariate gaussian copula with poisson marginals,

```
library(MASS)
fBivariateGaussianCopulaPoisson <- function(aVar,aCor,aMean,bMean){
    aSigma <- matrix(c(aVar,aCor,aCor,aVar),ncol = 2)
  vX <- mvrnorm(1,c(0,0),aSigma)
  vU <- pnorm(vX)
  aY <- qpois(vU[1],aMean)
  bY <- qpois(vU[2],bMean)
  return(c(aY,bY))
}

N <- 100
mCopula <- matrix(nrow = N,ncol=2)
for(i in 1:N){
    mCopula[i,] <- fBivariateGaussianCopulaPoisson(1,-0.5,5,3)  
}

```

I then use the below Stan file to attempt to estimate both the covariance matrix and the means of the marginals. It is able to reproduce the means of the marginals, but it does not seem to get the covariance matrix right. I have also tried the below but replacing inv\_Phi(mv[i,j]) with inv\_Phi(u[i,j]) but this does not seem to help. Note it is hard coded for the bivariate case.

Any ideas anyone?

```
functions{
  vector [] fFindBounds(row_vector vU, matrix L, int[] vY, vector vLambda){
    matrix[2,2] mBounds;
    vector[2] lLogProbU;
    matrix[2,2] u_star;
    vector[2] z_star;
    real z;
    vector[2] v;
    vector[2] out[2];

    mBounds[1,1] = poisson_cdf(vY[1]-1,vLambda[1]);
    mBounds[1,2] = poisson_cdf(vY[1],vLambda[1]);
    mBounds[2,1] = poisson_cdf(vY[2]-1,vLambda[2]);
    mBounds[2,2] = poisson_cdf(vY[2],vLambda[2]);

    z_star[1] = mBounds[1,1]/L[1,1];
    z_star[2] = mBounds[1,2]/L[1,1];

    u_star[1,1] = Phi(z_star[1]); 
    u_star[1,2] = Phi(z_star[2]);

    v[1] = u_star[1,1] + (u_star[1,2]-u_star[1,1]) * vU[1];
    z = inv_Phi(v[1]);

    u_star[2,1] = Phi((mBounds[2,1]-L[2,1] * z)/L[2,2]); 
    u_star[2,2] = Phi((mBounds[2,2]-L[2,1] * z)/L[2,2]); 

    v[2] = u_star[2,1] + (u_star[2,2]-u_star[2,1]) * vU[2];

    lLogProbU[1] = u_star[1,2] - u_star[1,1];
    lLogProbU[2] = u_star[2,2] - u_star[2,1];

    out[1] = lLogProbU;
    out[2] = v;

    return(out);
   }
}

data{
  int N;
  int Y[N,2];
}

transformed data{
  vector[2] vZeros;
  for(i in 1:2)
    vZeros[i] = 1;
}

parameters{
  cholesky_factor_cov[2] L; // covariance between Xs
  vector<lower=0>[2] lambda; // mean of each marginal poisson
  
  // augmented liKelihood parameters
  matrix<lower=0,upper=1>[N,2] u;
}

transformed parameters{
  matrix[N,2] x;
  vector[N] lLogProb;
  matrix[N,2] mLogProbU;
  matrix[N,2] z;
  vector[2] lOut[N,2];
  matrix[N,2] mv;
  
  for(i in 1:N){
    lOut[i] = fFindBounds(u[i],L,Y[i],lambda);
    mv[i] = to_row_vector(lOut[i,2]);
    mLogProbU[i] = to_row_vector(lOut[i,1]);
  }
    
  for(i in 1:N){
    for(j in 1:2){
      z[i,j] = inv_Phi(mv[i,j]);
    }
  }
      
      
  for(i in 1:N){
       x[i] = to_row_vector(L * to_vector(z[i]));
  }
    
  for(i in 1:N){
    lLogProb[i] = multi_normal_cholesky_lpdf(to_vector(x[i])|vZeros,L);
  }
}

model{
  for(i in 1:N){
    target += lLogProb[i];
    for(j in 1:2)
      target += log(mLogProbU[i,j]);
  }
  
  L ~ lkj_corr_cholesky(2.0);
  lambda ~ normal(0,10);
}

generated quantities{
  cov_matrix[2] Omega;
  Omega = L * L';
}
```

---

<div class="post-metadata">

**Author:** ![bgoodri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bgoodri/32/4451_2.png) [@bgoodri](https://discourse.mc-stan.org/u/bgoodri)\
**Post date:** [May 31, 2017, 4:13am UTC](https://discourse.mc-stan.org/t/gaussian-copula-for-discrete-marginals/691/10 "2017-05-31T04:13:14Z")

</div>

> [@Ben\_Lambert](#):
>
> cholesky\_factor\_cov[2] L; // covariance between Xs

That should be `cholesky_factor_corr[2]` for starters. Also, I don’t think `fFindBounds` is doing the right thing. When it calculates `zstar` I believe those should be values for standard normal variates, such that if you computed their standard normal CDF, those are the breakpoints in the quantile function for the Poisson. Specifically, the calculation of `z_star` should involve `inv_Phi`

So, if an outcome is 5 and its mean is 3, then

```no-highlight
z_star_high <- qnorm(ppois(4, 3))
qpois(pnorm(z_star_high), 3) # yields 4
qpois(pnorm(z_star_high) + 1e-13, 3) # now yields 5

z_star_low <- qnorm(ppois(4 - 1, 3))
qpois(pnorm(z_star_low), 3) # yields 3
qpois(pnorm(z_star_low) + 1e-13, 3) # now yields 4

```

So, you want a uniform between `u_star_low = Phi(z_star_low)` and `u_star_high = Phi(z_star_high)` but the Stan language only lets you put scalar bounds so you have a uniform between 0 and 1. So, do a change of variables there.

Then, for the second observed count, the critical values depend on `L` too.

---

<div class="post-metadata">

**Author:** ![Andre\_Pfeuffer](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/andre_pfeuffer/32/4939_2.png) [@Andre\_Pfeuffer](https://discourse.mc-stan.org/u/Andre_Pfeuffer)\
**Post date:** [May 31, 2017, 1:18pm UTC](https://discourse.mc-stan.org/t/gaussian-copula-for-discrete-marginals/691/11 "2017-05-31T13:18:20Z")

</div>

> ```
> real fcop(real tau,real mu,real ni) {
> return( pow( pow(mu,-tau) + pow(ni,-tau) - 1.0,-1/tau));
> 
> ```
> 
> }
> 
> real clayton\_poisson(int y1, int y2, real la1, real la2, real tau) {  
> if(y1==0 && y2==0) {  
> return fcop(tau,exp(-la1), exp(-la2));  
> }  
> else if(y1\>0 && y2==0) {  
> real C11m = poisson\_cdf(y1-1,la1);  
> real C1 = C11m + exp(poisson\_lpmf(y1 | la1));  
> return fcop(tau,C1, exp(-la2))  
> - fcop(tau,C11m, exp(-la2));  
> }  
> else if(y1==0 && y2\>0) {  
> real C21m = poisson\_cdf(y2-1, la2);  
> real C2 = C21m + exp(poisson\_lpmf(y2 | la2));  
> return fcop(tau,exp(-la1), C2)  
> - fcop(tau,exp(-la1), C21m);  
> }  
> else {  
> real C11m = poisson\_cdf(y1-1,la1);  
> real C1 = C11m + exp(poisson\_lpmf(y1 | la1));  
> real C21m = poisson\_cdf(y2-1, la2);  
> real C2 = C21m + exp(poisson\_lpmf(y2 | la2));  
> return fcop(tau, C1, C2)  
> - fcop(tau, C11m, C2)  
> - fcop(tau, C1, C21m)  
> + fcop(tau, C11m, C21m);  
> }  
> }

la1, la2 is on log-scale. It uses a bivariate poisson using clayton copula. The gaussian copula I posted one time in the  
google Stan discussion board.

Added: Take a look at A PRIMER ON COPULAS FOR COUNT DATA, CHRISTIAN GENEST  
AND JOHANNA NESLEHOVA for details and how and if you can expand this in the  
multivariate case. There is also a multivariate poisson distribution.

Andre

---

<div class="post-metadata">

**Author:** ![Ben\_Lambert](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ben_lambert/32/4407_2.png) [@Ben\_Lambert](https://discourse.mc-stan.org/u/Ben_Lambert)\
**Post date:** [June 2, 2017, 1:29am UTC](https://discourse.mc-stan.org/t/gaussian-copula-for-discrete-marginals/691/12 "2017-06-02T01:29:37Z")

</div>

Hi Ben,

Thanks for your advice here.

I agree I need to be taking the inverse standard normal of the poisson cdf value first. I have also changed to a correlation cholesky matrix.

Now I have the following code, which unfortunately seems to generate ```Rejecting initial value: gradient evaluated at the initial value is not finite’. I have been debugging by print but can’t seem to find the cause of this. Any ideas?

```
functions{
  vector [] fFindBounds(row_vector vU, matrix L, int[] vY, vector vLambda){
    matrix[2,2] mBounds;
    vector[2] lLogProbU;
    vector[2] u_star_lower;
    vector[2] u_star_upper;
    real z_star_lower;
    real z_star_upper;
    real z;
    vector[2] v;
    vector[2] out[2];

mBounds[1,1] = inv_Phi(poisson_cdf(vY[1]-1,vLambda[1]));
mBounds[1,2] = inv_Phi(poisson_cdf(vY[1],vLambda[1]));
mBounds[2,1] = inv_Phi(poisson_cdf(vY[2]-1,vLambda[2]));
mBounds[2,2] = inv_Phi(poisson_cdf(vY[2],vLambda[2]));

z_star_lower = mBounds[1,1]/L[1,1];
z_star_upper = mBounds[1,2]/L[1,1];

u_star_lower[1] = Phi(z_star_lower); 
u_star_upper[1] = Phi(z_star_upper);

v[1] = u_star_lower[1] + (u_star_upper[1]-u_star_lower[1]) * vU[1];

z = inv_Phi(v[1]);

u_star_lower[2] = Phi((mBounds[2,1]-(L[2,1] * z))/L[2,2]); 
u_star_upper[2] = Phi((mBounds[2,2]-(L[2,1] * z))/L[2,2]); 

v[2] = u_star_lower[2] + (u_star_upper[2]-u_star_lower[2]) * vU[2];

lLogProbU[1] = u_star_upper[1] - u_star_lower[1];
lLogProbU[2] = u_star_upper[2] - u_star_lower[2];

out[1] = lLogProbU;
out[2] = v;

return(out);
  }
}

data{
  int N;
  int Y[N,2];
}

transformed data{
  vector[2] vZeros;
  for(i in 1:2)
    vZeros[i] = 0;
}

parameters{
  cholesky_factor_corr[2] L; // correlation between Xs
  vector<lower=0>[2] lambda; // mean of each marginal poisson
  
  // augmented liKelihood parameters
  matrix<lower=0,upper=1>[N,2] u;
}

transformed parameters{
  matrix[N,2] x;
  vector[N] lLogProb;
  matrix[N,2] mLogProbU;
  matrix[N,2] z;
  vector[2] lOut[N,2];
  matrix[N,2] mv;
  
  for(i in 1:N){
    lOut[i] = fFindBounds(u[i],L,Y[i],lambda);
    mv[i] = to_row_vector(lOut[i,2]);
    mLogProbU[i] = to_row_vector(lOut[i,1]);
  }
    
  for(i in 1:N){
    for(j in 1:2){
      z[i,j] = inv_Phi(mv[i,j]);
    }
  }
      
  for(i in 1:N){
       x[i] = to_row_vector(L* to_vector(z[i]));
  }
    
  for(i in 1:N){
    lLogProb[i] = multi_normal_cholesky_lpdf(to_vector(x[i])|vZeros,L);
  }
}

model{
  for(i in 1:N){
    target += lLogProb[i];
    for(j in 1:2)
      target += log(mLogProbU[i,j]);
  }
  
  L ~ lkj_corr_cholesky(2.0);
  lambda ~ normal(0,10);
}
```

---

<div class="post-metadata">

**Author:** ![Ben\_Lambert](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ben_lambert/32/4407_2.png) [@Ben\_Lambert](https://discourse.mc-stan.org/u/Ben_Lambert)\
**Post date:** [June 5, 2017, 4:36pm UTC](https://discourse.mc-stan.org/t/gaussian-copula-for-discrete-marginals/691/13 "2017-06-05T16:36:34Z")

</div>

Hi Ben,

Sorry to bother you, but am really quite stuck on this one. Any ideas on the cause of the issue I mentioned in the previous message?

Best,

Ben

---

<div class="post-metadata">

**Author:** ![bgoodri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bgoodri/32/4451_2.png) [@bgoodri](https://discourse.mc-stan.org/u/bgoodri)\
**Post date:** [June 6, 2017, 3:57pm UTC](https://discourse.mc-stan.org/t/gaussian-copula-for-discrete-marginals/691/14 "2017-06-06T15:57:40Z")

</div>

On Discourse, you have to tag @bgoodri in order for it to show up in the list at the top right. Also, for things like this, you have to install StanHeaders and RStan from the 2.14.x release in order to see the actual error message, which in your case is

```
[407] "Exception thrown at line 88: multi_normal_cholesky_log: Random variable[1] is -nan, but must not be nan!"  

```

It looks like `inv_Phi` is returning NaN, which I would guess means sometimes its input is outside of the (0,1) interval.

---

<div class="post-metadata">

**Author:** ![Ben\_Lambert](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ben_lambert/32/4407_2.png) [@Ben\_Lambert](https://discourse.mc-stan.org/u/Ben_Lambert)\
**Post date:** [June 7, 2017, 12:00am UTC](https://discourse.mc-stan.org/t/gaussian-copula-for-discrete-marginals/691/15 "2017-06-07T00:00:59Z")

</div>

@bgoodri – thanks for your response. I can replicate the error for the few datasets for which the code appears to run. These are inevitably the smaller (fewer N) datasets, although it only appears to run about 1 time in 10, for a sample size of N=5. It does appear to produce reasonable results for the times when it does execute. However I cannot figure out what drives the failure during the 9/10 times when it fails, apart from it seems more likely to run when N is small.

There are occasions when the inputs appear to print as either 0 or 1 (for example when printing v[1]), and hence inv\_Phi evaluates to - or + inf. However I cannot understand what is driving the sampler to visit these values…

I have tried passing initial values for the parameters which are sensible but this doesn’t seem to ameliorate the situation either.

Am really at wit’s end here, so any help would be most appreciated.

Best,

Ben

---

<div class="post-metadata">

**Author:** ![bgoodri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bgoodri/32/4451_2.png) [@bgoodri](https://discourse.mc-stan.org/u/bgoodri)\
**Post date:** [June 7, 2017, 2:45am UTC](https://discourse.mc-stan.org/t/gaussian-copula-for-discrete-marginals/691/16 "2017-06-07T02:45:15Z")

</div>

The observed counts that are zero mess up `inv_Phi(poisson_cdf(vY[1]-1,vLambda[1]))` and `inv_Phi(poisson_cdf(vY[2]-1,vLambda[2]))`.  
You need to handle those cases separately because it is just uniform between 0 and whatever the upper bound is.

---

<div class="post-metadata">

**Author:** ![Ben\_Lambert](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ben_lambert/32/4407_2.png) [@Ben\_Lambert](https://discourse.mc-stan.org/u/Ben_Lambert)\
**Post date:** [June 8, 2017, 6:58pm UTC](https://discourse.mc-stan.org/t/gaussian-copula-for-discrete-marginals/691/17 "2017-06-08T18:58:12Z")

</div>

@bgoodri Thank you very much for all your help with this. I have now coded up the example for a multivariate distribution with K dimensions (see below). All seems to work – the model is actually excellent at returning back the parameters that I use to generate the data (R code also below). Many thanks again for your help here, and include the code below in case others find it useful. I am, in particular, looking forward to extending this model to encompass a variety of discrete marginals: bernoulli, binomial, and negative binomial for starters. Best, Ben

R code to generate data from a Gaussian Copula with Poisson marginals (using a covariance matrix aSigma, and with marginal means given by the vector lMeans):

```
fMultivariateGaussianCopulaPoisson <- function(aSigma,lMeans){
  K <- length(lMeans)
  vX <- mvrnorm(1,rep(0,K),aSigma)
  vU <- pnorm(vX)
  lY <- vector(length=K)
  for(i in 1:K)
    lY[i] <- qpois(vU[i],lMeans[i])
  return(lY)
}

// For N observations from the above
fGenerateMultivariateCopulaPoisson <- function(N,aSigma,lMeans){
  K <- nrow(aSigma)
  mCopula <- matrix(nrow = N,ncol=K)
  for(i in 1:N){
    mCopula[i,] <- fMultivariateGaussianCopulaPoisson(aSigma,lMeans) 
  }
  return(mCopula)
}

```

Stan code (note I had to constrain the Poisson marginal parameters to avoid initialising these parameters with too large values that caused the log probability to go to -inf):

```
functions{
  vector [] fFindBounds(row_vector vU, matrix L, int[] vY, vector vLambda){
    int K = cols(vU);
    matrix[K,2] mBounds;
    vector[K] lLogProbU;
    vector[K] u_star_lower;
    vector[K] u_star_upper;
    real z_star_lower;
    real z_star_upper;
    vector[K] z;
    vector[K] v;
    vector[K] out[2];
    for(k in 1:K){
      int km1 = k - 1;
      if(vY[k]<0.001){ ## upper bound only
        mBounds[k,1] = -1;
        mBounds[k,2] = inv_Phi(poisson_cdf(vY[k],vLambda[k]));
        z_star_lower = -1;
        z_star_upper = (mBounds[k,2] - ((k > 1) ? L[k,1:km1] * head(z,km1) : 0)) / L[k,k];
        u_star_lower[k] = 0;
        u_star_upper[k] = Phi(z_star_upper);
      }else{
        mBounds[k,1] = inv_Phi(poisson_cdf(vY[k]-1,vLambda[k]));
        mBounds[k,2] = inv_Phi(poisson_cdf(vY[k],vLambda[k]));
        z_star_lower = (mBounds[k,1] - ((k > 1) ? L[k,1:km1] * head(z,km1) : 0)) / L[k,k];
        z_star_upper = (mBounds[k,2] - ((k > 1) ? L[k,1:km1] * head(z,km1) : 0)) / L[k,k];
        u_star_lower[k] = Phi(z_star_lower); 
        u_star_upper[k] = Phi(z_star_upper);
      }
      v[k] = u_star_lower[k] + (u_star_upper[k]-u_star_lower[k]) * vU[k];
      z[k] = inv_Phi(v[k]);
      lLogProbU[k] = u_star_upper[k] - u_star_lower[k];
    }
    
    out[1] = lLogProbU;
    out[2] = v;
    
    return(out);
  }
}

data{
  int N;
  int K;
  int Y[N,K];
}

transformed data{
  vector[K] vZeros;
  for(i in 1:K)
    vZeros[i] = 0;
}

parameters{
  cholesky_factor_corr[K] L; // correlation between Xs
  vector<lower=0,upper=30>[K] lambda; // mean of each marginal poisson
  vector<lower=0>[K] sigma;
  
  // augmented likelihood parameters
  matrix<lower=0,upper=1>[N,K] u;
}

transformed parameters{
  matrix[N,K] x;
  vector[N] lLogProb;
  matrix[N,K] mLogProbU;
  matrix[N,K] z;
  vector[K] lOut[N,2];
  matrix[N,K] mv;
  cholesky_factor_cov[K] L_cov;
  
  L_cov = diag_pre_multiply(sigma,L);
  
  for(i in 1:N){
    lOut[i] = fFindBounds(u[i],L_cov,Y[i],lambda);
    mv[i] = to_row_vector(lOut[i,2]);
    mLogProbU[i] = to_row_vector(lOut[i,1]);
  }
    
  for(i in 1:N){
    for(j in 1:K){
      z[i,j] = inv_Phi(mv[i,j]);
    }
  }
      
  for(i in 1:N){
       x[i] = to_row_vector(L_cov * to_vector(z[i]));
  }
    
  for(i in 1:N){
    lLogProb[i] = multi_normal_cholesky_lpdf(to_vector(x[i])|vZeros,L_cov);
  }
}

model{
  for(i in 1:N){
    target += lLogProb[i];
    for(j in 1:K)
      target += log(mLogProbU[i,j]);
  }
  
  L ~ lkj_corr_cholesky(2.0);
  lambda ~ normal(5,10);
  sigma ~ lognormal(0,1);
}

generated quantities{
  cov_matrix[K] Omega;
  Omega = L_cov * L_cov';
}
```

---

<div class="post-metadata">

**Author:** ![bgoodri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bgoodri/32/4451_2.png) [@bgoodri](https://discourse.mc-stan.org/u/bgoodri)\
**Post date:** [June 8, 2017, 7:56pm UTC](https://discourse.mc-stan.org/t/gaussian-copula-for-discrete-marginals/691/18 "2017-06-08T19:56:44Z")

</div>

> [@Ben\_Lambert](#):
>
> Stan code (note I had to constrain the Poisson marginal parameters to avoid initialising these parameters with too large values that caused the log probability to go to -inf):

Did it not work to just initialize the parameters with smaller values?

---

<div class="post-metadata">

**Author:** ![Ben\_Lambert](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ben_lambert/32/4407_2.png) [@Ben\_Lambert](https://discourse.mc-stan.org/u/Ben_Lambert)\
**Post date:** [June 8, 2017, 7:58pm UTC](https://discourse.mc-stan.org/t/gaussian-copula-for-discrete-marginals/691/19 "2017-06-08T19:58:23Z")

</div>

Yes, this is what I actually have done subsequently.

---

<div class="post-metadata">

**Author:** ![anhsmith](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/anhsmith/32/10378_2.png) [@anhsmith](https://discourse.mc-stan.org/u/anhsmith)\
**Post date:** [February 19, 2019, 2:25am UTC](https://discourse.mc-stan.org/t/gaussian-copula-for-discrete-marginals/691/20 "2019-02-19T02:25:51Z")

</div>

@Ben_Lambert thank you for posting your code here. I’ve tried it out and it works well. I work a lot with counts of organisms and these sorts of models could be extremely useful. I’d be interested to hear if you’ve developed this any further, e.g., for other marginal distributions such as the negative binomial, as you suggested.

[Next page](https://discourse.mc-stan.org/t/gaussian-copula-for-discrete-marginals/691.md?page=2)