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';
}
```