How simulate correlated Poisson distributions

12 次查看（过去 30 天）
Pete sherer
Pete sherer 2018-10-7
评论： 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.

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

采纳的回答

Bruno Luong
Bruno Luong 2018-10-7
编辑：Bruno Luong 2018-10-8
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更早的评论
Pete sherer
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
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
Jeff Miller 2018-10-8
编辑：Jeff Miller 2018-10-8
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更早的评论
Jeff Miller
Jeff Miller 2018-10-15
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
Jeff Miller 2018-10-16
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
Pete sherer 2018-10-17
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
Pete sherer 2018-10-17
also
nb = NegativeBinomial(alpha,p);
Undefined function or variable 'NegativeBinomial'.
Jeff Miller
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Binomial Distribution 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by