5 views (last 30 days)

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.

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.

Opportunities for recent engineering grads.

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

Start Hunting!
## 0 Comments

Sign in to comment.