How simulate correlated Poisson distributions
17 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Pete sherer
am 7 Okt. 2018
Kommentiert: Jeff Miller
am 17 Okt. 2018
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 Kommentare
Akzeptierte Antwort
Bruno Luong
am 7 Okt. 2018
Bearbeitet: Bruno Luong
am 8 Okt. 2018
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 Kommentare
Bruno Luong
am 14 Okt. 2018
Bearbeitet: Bruno Luong
am 14 Okt. 2018
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.
Weitere Antworten (2)
Jeff Miller
am 8 Okt. 2018
Bearbeitet: Jeff Miller
am 8 Okt. 2018
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 Kommentare
Jeff Miller
am 15 Okt. 2018
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
am 16 Okt. 2018
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
am 17 Okt. 2018
2 Kommentare
Jeff Miller
am 17 Okt. 2018
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.
Siehe auch
Kategorien
Mehr zu Binomial Distribution finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!