How simulate correlated Poisson distributions

Hi Is there a way to simulate correlated RVs where each RV follows poisson distribution? I have 2 RVs X1 and X2, and both follow Poisson distribution. I would like simulate final results such that I can control correlation between X1 and X2. Either simulating straight correlated process or post-processing of 2 independent run would work.

Answer 1

[Bruno Luong]

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)

[Bruno Luong]

I haven't though how to express the correlation of the induced Poisson RV from their induced uniform distributions. There might be an explicit formula.

[Pete sherer]

Hi Bruno, Is this typical to see the desired Xcorr of 0.6, but the end correlation is reduced to just under 0.4?

[Bruno Luong]

First it's closer to 0.5 than 0.4.

Second: As I said earlier in my subsequent comment and in the code (look at the comment above Xcorr set instruction), the Xcorr is not the correlation of the final Poisson RV but the correlation of the uniform distributions used to generated the Poisson. They are related monotonically but I haven't workout the explicit formula.

Still it gives you a way to control the correlation.

[Bruno Luong]

I couldn't find a rigorous close formula related the correlation between Y and X. So I derive a completely empirical formula:

n = 100000;
lambda  = 10; % scalar or vector of 1x2

Ycorr = 0.6;
% here is the emperical formula
alpha = polyval([0.0708 0.6408], log(lambda));
Xcorr = 1 - (1-Ycorr).^(1/alpha);

M12 = 2*sin(pi*Xcorr/6);
M = [1, M12;
    M12, 1];
C = chol(M);

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:3*lambda])
corrcoef(Y)
mean(Y)

[Bruno Luong]

The code could be accelerate by generating a bunch of X before the while loop. The case lambda >> 1 (let's say > 10) can also be better handle, using Gaussian approximation etc... this is well documented in textbook.

[Pete sherer]

Thank you very much for providing this insightful script. If possible, could you comments on?

• To make this simulation more generic and not limited to correlated poisson. In theory I want to be able to mix poisson with poisson-mixture (e.g., Poisson Gamma). I could generate those poisson and poisson-gamma indepedently beforehands.
• For more than 2 RVs, the C should be C= C/sqrt( nRVs), right?

X % [n, nRVs] presimulated independent RVs   
C = C / sqrt( nRVs);
corrX= X* C;

[Bruno Luong]

The factor sqrt(2) comes from the relation

    normcdf(x)=1/2*erfc(−x/sqrt(2)).

So you shouldn't change it.

For extension to Poisson Gamma; I guess you can generate N-lambda independently from alpha/beta, and make the correlation works from X alone.

[Bruno Luong]

Here is the code for more than 2 RVs and might be a way for Gamma-Poisson. But after playing with it, depending on parameters one might need to generate a correlated gamma-distribution as well. The on/off approach is too violent and it can create too correlated or not enough correlated outcome.

% number of RVs
nvars = 3;

% number of random samples per variables
n = 100000;

% Desired correlation
ycorrtab = [0.3;  % y12
            0.2;  % y13
            0.1]; % y23

% Fix lambda==true -> Poisson distribution with lambda = alpha*beta
% Otherwise Gamma-Poisson with parameters alpha and beta
fixlambda = true;
alpha = 8;
beta = 0.5;
% If TRUE we use the same lambda for all RV
% they will be stronger correlated
commonlambda = true; %all(ycorrtab > 1/3);

% Allocate array
Y = zeros(n,nvars);
lambda = alpha*beta;
Ycorr = ycorrtab(:).';

if fixlambda
    % empirical correlation transformation
    gamma = polyval([0.0708 0.6408], log(lambda));
    Xcorr = 1 - (1-Ycorr).^(1./gamma);
else
    % WARNING: Conversion not valid for gamma-Poisson
    Xcorr = 1*Ycorr;
end

xcorr = 2*sin(pi*Xcorr/6);
M = zeros(nvars);
M(triu(true(nvars),1)) = xcorr;
M = M+M';
M(1:nvars+1:end) = 1;
C = chol(M);
C = C / sqrt(2);

for r=1:n
    if ~fixlambda
        if commonlambda
            % same lambda for all RV: stronger correlation
            lambda = gamrnd(alpha,beta,1,1);           
        else
            % different lambda for all VR: weaker correlation
            lambda = gamrnd(alpha,beta,1,nvars);
        end
    end
    % Knuth algorithm
    k = zeros(1,nvars);
    p = ones(1,nvars);
    L = exp(-lambda);    
    b = p > L;
    while any(p > L)
        k = k+double(b);
        X = (erf(randn(1,nvars)*C) + 1) / 2;
        p = p .* X;
        b = p > L;
    end
    Y(r,:) = k-1;
end

close all
hist(Y,[0:4*alpha*beta])
c = corrcoef(Y)
mean(Y)
if fixlambda
    title(sprintf('Poisson \\lambda=%g', lambda));
else
    title(sprintf('\\Gamma-Poisson \\alpha=%g, \\beta=%g', alpha, beta));
end

[Pete sherer]

Thanks very much Bruno. This is greatly insightful.

[Pete sherer]

Hi Bruno, By the way I was trying your code with alpha(shape)=1.0 and beta=2.0. I was hoping that i will get regular poisson back, but it seems not the case when comparing output with poissrnd.

[Bruno Luong]

No, I use the definition as I stated above, to get the poisson one must have (

• (beta -> 0)
• (alpha*beta = lambda), so alpha -> infinity

In your test, I would put beta=0.1, alpha=20, that would simulates RV similar to Poisson with lambda=2.

Answer 2

[Jeff Miller]

You can do this with the RandGen class in Cupid . The code will look something like this:

NRVs = 2;
RVs = cell(NRVs,1);
RVs{1}=Poisson(23);   % The Poisson parameter is the mean
RVs{2}=Poisson(12);
TargetCorrs = [1  0.3; 0.3 1];  % Replace 0.3 with whatever correlation you want.
mymultivar = RandGen(RVs,TargetCorrs,'Adjust','NSteps',500);
r=mymultivar.Rands(10000);  % r is the array of random numbers with desired marginals & correlation

Look at DemoRandGen.m for more details.

[Pete sherer]

Hi Jeff, I ran your example and got this missing function error msg.

Undefined function or variable 'ExtractNamei'.

[Torsten]

https://github.com/milleratotago/RawRT/blob/master/ExtractNamei.m

[Jeff Miller]

Yes, Cupid uses some routines from ExtractNameVal, so you also need that. I guess I need to document that more prominently--sorry.

[Pete sherer]

Thanks very much Jeff.

[Pete sherer]

Undefined function or variable 'LegalCorrMatrix'.

Error in RandGen/Init (line 75) assert(LegalCorrMatrix(obj.RhoControllers),'Error! The requested correlation matrix is impossible.!');

Error in RandGen (line 57) Init(obj,varargin{1},varargin{2},varargin{3:end});

[Jeff Miller]

Sorry again. Here it is:

function [ bool ] = LegalCorrMatrix( CorrMatrix )
% Returns 1 if the input matrix is a legal correlation matrix,
%         0 otherwise.

[~, err] = cholcov(CorrMatrix);
if err == 0
    bool = true;
else
    bool = false;
end

end

[Pete sherer]

Thanks. I can run through now. Does your code support Poisson-Gamma simulation?

[Jeff Miller]

Good that it is working.

I'm not entirely sure what you mean by "Poisson-Gamma simulation". According to Wikipedia, it looks like you would just need the negative binomial distribution for that. Is that what you are after, or would you need something more?

[Bruno Luong]

As I understood, Poisson-Gamma is the probability law of the number event occurring in a unit time, but the average rate (lambda) of the events instead fix as with Poisson, follows a Gamma law with alpha/beta parameters.

[Pete sherer]

I was thinking of Poisson-Mixture in general. Don't worry, I can generate the marginal with that.

[Jeff Miller]

@Bruno: Yes, that was my understanding too. But according to Wikipedia, the probability law resulting from that process seems to be a negative binomial.

@Pete: Cupid can make random variables that are any arbitrary mixtures of (any number of) Poissons, e.g.:

rvmix = Mixture(0.3,Poisson(4),0.35,Poisson(5),0.35,Poisson(7.2))

and you can then use one of these distributions as one of your marginals. It might be a bit slow, though.

[Bruno Luong]

Thanks I'll read wikipedia. ;-)

Poisson Mixture would be indeed easier to deal with than jumping around as Gamma.

[Pete sherer]

@Jeff, I was naive thinking that I can replace your RVs{} with pre-simulated Pois-Gamma. However it seems your RVs is an object that specified distribution. Is there a way I can hijack that step and only passing pre-simulated data before your code doing the correlation among RVs? Or is there a way to convert function below to your RV object?

function simVec= PoisGamma( alpha, beta, nSim)
simLambda = gamrnd( alpha, beta, [nSim, 1]);
simVec= poissrnd( simLambda, [nSim, 1]);

[Jeff Miller]

No, you can't pass data, but I think you can get what you want anyway using the NegativeBinomial distribution (which I added to Cupid yesterday). Again, I refer you to Wikipedia which explains the equivalence. Here is a small script illustrating the equivalence:

alpha=4;
beta=22;
nSim=100000;
simLambda = gamrnd( alpha, beta, [nSim, 1]);
simVec= poissrnd( simLambda, [nSim, 1]);
figure; histogram(simVec)

% Equivalent negative binomial, according to https://en.wikipedia.org/wiki/Negative_binomial_distribution
p = 1/(1+beta);
nb = NegativeBinomial(alpha,p);
simVec2 = nb.Random([nSim,1]);
figure; histogram(simVec2);

For all values of alpha and beta that I have tried, the histograms look identical, so this might convince you of the equivalence even if the math in Wikipedia does not. I think this means you can get what you want using the NegativeBinomial, which is a regular cupid distribution.

[Jeff Miller]

Another option--just added--is to use a "List" random variable, which is defined by an arbitrary vector of discrete random values and another arbitrary vector of the probabilities of each value. You would get these values by random generation, more or less like this:

simLambda = gamrnd( alpha, beta, [nSim, 1]);
simVec= poissrnd( simLambda, [nSim, 1]);
[Xs, ~, addresses] = unique(simVec);
Counts = accumarray(addresses,1);
Ps = Counts / sum(Counts);
gpRV = List(Xs',Ps');  % Note transpose relative to unique & accumarray outputs

Now gpRV is a Cupid distribution object that you can assign to RV{1}, etc, and use with RandGen. It won't give you exactly the values in simVec, but it will give you values from the marginal distribution defined by simVec.

Answer 3

[Pete sherer]

Is the List() function a default one? I got this error message

gpRV = List(Xs',Ps'); % Note tran
Undefined function 'List' for input arguments of type 'double'.

[Pete sherer]

also

nb = NegativeBinomial(alpha,p);
Undefined function or variable 'NegativeBinomial'.

[Jeff Miller]

List and NegativeBinomial were both added to the repo since we started this discussion. If you just downloaded Cupid once at the beginning, you didn't get them, and you'll need to update to the latest version for these functions to be defined.