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

Question

Daniel Niu 2022-10-18

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1829703-what-is-the-problem-of-my-code-doing-inverse-poisson-cumulative-distribution

评论： Daniel Niu 2022-10-19

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

2 个评论
显示无隐藏无

Torsten 2022-10-18

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 2022-10-18

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 2022-10-18

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1829703-what-is-the-problem-of-my-code-doing-inverse-poisson-cumulative-distribution#answer_1078193

编辑：Torsten 2022-10-18

在 MATLAB Online 中打开

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

3 个评论
显示 1更早的评论隐藏 1更早的评论

Torsten 2022-10-19

在 MATLAB Online 中打开

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 2022-10-19

Thank you very much!

请先登录，再进行评论。

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

2 个评论
显示无隐藏无

采纳的回答

3 个评论
显示 1更早的评论隐藏 1更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

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

2 个评论 显示 无隐藏 无

采纳的回答

3 个评论 显示 1更早的评论隐藏 1更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

2 个评论
显示无隐藏无

3 个评论
显示 1更早的评论隐藏 1更早的评论