# Generating random draws from a specific poisson distribution

5 views (last 30 days)
ektor on 23 Oct 2015
Edited: ektor on 25 Oct 2015
Dear all,
I have the following variation of the poisson distribution
P(Y=y)=((mean/vbn)^(y))*(nv^(y-1))*(exp(-mean*nv/vbn))/factorial(y);
where y is the count, vbn=1+overdisp*mean where overdisp is a dispersion parameter that can be < or >0 and nv= 1+overdisp*y;
I want to generate random draws from this distribution. So I wrote the following code
function yy = Poissaw(mean,ns, overdisp);
YY= zeros(ns,1);
qqq=(0:100)';
vbn=1+overdisp*mean;
nv= 1+overdisp.*qqq
for i = 1:ns
fff= ((mean/vbn).^(qqq)).*(nv.^(qqq-1)).*(exp(-mean.*nv./vbn))./factorial(qqq);
ff=cumsum(fff);
YY(i)= sum( rand > ff);
end
Could some verify the correctness (or not) of the above function?
Many thanks

Walter Roberson on 24 Oct 2015
The body of your "for" loop does not involve "i" except as the index to store the value. The variable "fff" that you are calculating does not involve "i" and there are no changes in the loop to any variable that you are calculating from. Therefore the value of "fff" and "ff" that you calculate each loop is the same. Therefore you can calculate them outside the loop.
Inside the loop you are essentially doing a "find" operation, counting the number of items before the random value is le ff. You can vectorize that.
YY = sum(bsxfun(@gt, rand(ns,1), ff.'));
It is not immediately obvious to me that the final value in ff is 1.0; it is not obvious to me that it could not be greater than or less than 1.0 . It is therefor not clear that rand() is proper here. If it is, then if the distribution has an infinite tail, you appear to be truncating the distribution at 100. It is not obvious that is justified numerically.
If you do decide to expand to a wider range for qqq, then you will have the ultimate limit that rand() cannot produce any value between (1-2^(-53)) and 1 because there simply are no representable double precision numbers in that range . On the lower end of the scale, most of the random number generators cannot produce a rand() value less than 2^(-53). There is, however, one random number generator that can produce any finite double precision number (including greater than 1), but you would need to select it specifically.

Show 1 older comment
ektor on 24 Oct 2015
Any help?
Walter Roberson on 24 Oct 2015
nv is not visible to the function. You do not have an "end" matching the "function" statement, so it cannot be the case that you are using a nested function.
Naming a variable "mean" is confusing and interferes with the use of "mean" as the MATLAB function.
ektor on 24 Oct 2015
Thank you very much Walter,
So now I have
function yy = Poissaw(meank,ns, overdisp,nmax);
qqq=(0:nmax)';
nv= 1+overdisp.*qqq
YY= zeros(ns,1);
for i = 1:ns
vbn=1+overdisp*meank(i);
fff= ((meank(i)/vbn).^(qqq)).*(nv.^(qqq-1)).*(exp(-meank(i).*nv./vbn))./factorial(qqq);
ff=cumsum(fff);
YY(i)= sum( rand > ff);
end
end
But even in this case, I still obtain ff values that do not sum to 1. Any suggestion/code on how to sample from this distribution? I am stuck and do not know how to move on.
Could you, please, provide some sample code?