Generating random draws from a specific poisson distribution

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

Answer 1

[Walter Roberson]

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.

[ektor]

Thanks Walter.

You are right that i does not involve inside the loop. Assuming that "mean" varies across "i" (so "mean" is a vector) the function becomes

 qqq=(0:100)';
nv= 1+overdisp.*qqq 
 function yy = Poissaw(mean,ns, overdisp);
YY= zeros(ns,1);
 for  i = 1:ns  
vbn=1+overdisp*mean(i);
fff= ((mean(i)/vbn).^(qqq)).*(nv.^(qqq-1)).*(exp(-mean(i).*nv./vbn))./factorial(qqq);
ff=cumsum(fff);
 YY(i)= sum( rand > ff);

 end

Would that be correct?

You are also correct in that that "ff" could be "inf" because I just noticed that ff takes such values. HOw could I make ff to be at maximum 1.?

Or more generally how would sample from this distribution?

Any help is greatly appreciated. Many thanks

[ektor]

Any help?

[Walter Roberson]

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]

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?