what is the problem of my code doing inverse Poisson cumulative distribution?

Dear friends,

I write a code doing inverse Poisson cumulative distribution. However, the code get same results compared to icdf function when u>0.5. When u<0.5, the result is 1 less than the result from icdf.

function m = iPoisson(lambda,u)

p(1)=exp(-lambda)/factorial(0);

absolute=zeros(1,1000);

summation=p(1);

%u=0.6322;

%u=rand(1);

for k=2:1000

p(k)=lambda.^(k-1)*exp(-lambda)./factorial(k-1);

summation=summation+p(k);

distance(k-1)=abs(summation-u);

[a,b]=min(distance);

m=b;

end

end

[Torsten]

Use a while-loop instead of a for loop for which you don't know when to stop it (do you know how large factorial(1000) is ?) .

[Daniel Niu]

Dear Torsten,

Thank you for your reply. I don't know how to use while because I don't konw how to define

k in while loop

Your help would be highly appreciated!

Best regards!

Answer 1

[Torsten]

lambda = 10;
u = 0.6;
m = iPoisson(lambda,u)

m = 11

m1 = icdf('Poisson',u,lambda)

m1 = 11

function m = iPoisson(lambda,u)
  if u == 1
    m = Inf;
    return
  end
  s = 1;
  m = 0;
  quot = 1;
  while s*exp(-lambda) < u 
    m = m + 1;
    quot = quot*lambda/m;
    s = s + quot;
  end
end

[Daniel Niu]

Dear Torsten,

The code works. Thank you very much!

However, I can not understand the mathematics. Why there is no factorial in the code? And how to get the m is the most right value just by s*exp(-lambda) < u?

May be s*exp(-lambda)>u a little, and is the nearest value

I can not figure out that.

Your help would be highly appreciated!

[Torsten]

You search for the first value m for which

exp(-lambda) * sum_{i=0}^{i=m} lambda^i/factorial(i) >= u. 

Since

s(k) = sum_{i=0}^{i=k} lambda^i/factorial(i)

is increasing with k, you need to add the terms

quot(i) = lambda^i / factorial(i)

by writing

s = s + quot

as long as

exp(-lambda)*s(k) < u.

If

exp(-lambda)*s(k) >= u

for the first time, then

m = k-1

Since

quot(i) = lambda/i * quot(i-1)

with

quot(0) = 1

you can update

quot = quot * lambda/i

with an initial value for quot

quot = 1

and you do not need to recalculate

 quot = lambda^i / factorial(i) 

for every value of i. Moreover, the simple update

 quot = quot * lambda/i

prevents overflow in numerator and/or denominator since lambda^i and factorial(i) both can become very large.

[Daniel Niu]

Thank you very much!