How simulate correlated Poisson distributions

20 次查看（过去 30 天）

Pete sherer 2018-10-7

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/422662-how-simulate-correlated-poisson-distributions

评论： Jeff Miller 2018-10-17

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.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Bruno Luong 2018-10-7

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/422662-how-simulate-correlated-poisson-distributions#answer_340200

编辑：Bruno Luong 2018-10-8

在 MATLAB Online 中打开

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)

11 个评论
显示 9更早的评论隐藏 9更早的评论

Bruno Luong 2018-10-9

编辑：Bruno Luong 2018-10-9

在 MATLAB Online 中打开

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 2018-10-14

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 2018-10-14

编辑：Bruno Luong 2018-10-14

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.

请先登录，再进行评论。

更多回答（2 个）

Jeff Miller 2018-10-8

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/422662-how-simulate-correlated-poisson-distributions#answer_340257

编辑：Jeff Miller 2018-10-8

在 MATLAB Online 中打开

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.

15 个评论
显示 13更早的评论隐藏 13更早的评论

Jeff Miller 2018-10-15

在 MATLAB Online 中打开

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 2018-10-16

在 MATLAB Online 中打开

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.

请先登录，再进行评论。

Pete sherer 2018-10-17

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/422662-how-simulate-correlated-poisson-distributions#answer_341868

在 MATLAB Online 中打开

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'.

2 个评论
显示无隐藏无

Pete sherer 2018-10-17

在 MATLAB Online 中打开

also

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

Jeff Miller 2018-10-17

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.

请先登录，再进行评论。

请先登录，再回答此问题。