How simulate correlated Poisson distributions

One way is to apply Knuth's Poisson generation using two correlated uniform distribution.

n = 10000;

% correlation coef of uniform distribution
% NOT of the poisson, but they are monotonically related
Xcorr = 0.6;
M12 = 2 * sin(pi * Xcorr / 6);
M = [1, M12;
     M12, 1];
C = chol(M);

% expectation value(s)
lambda  = 4; % scalar or vector of 1x2
Y = zeros(n,2);
L = exp(-lambda);
C = C / sqrt(2);
for r=1:n
    k = zeros(1,2);
    p = ones(1,2);
    b = p > L;
    while any(p > L)
        k = k+double(b);
        X = (erf(randn(1,2)*C) + 1) / 2;
        p = p .* X;
        b = p > L;
    end
    Y(r,:) = k-1;
end

hist(Y,[0:12])
corrcoef(Y)
mean(Y)