Poisson Function

[topkatz]

Hi.

For a given mean, m, I have computed individual Poisson probabilities three ways.

1. Using the SAS POISSON function. This is cumulative, so
P(m,n) = POISSON(m,n) - POISSON(m,(n-1))

2. Using the formula
P(m,n) = exp(-m) * (m**n) / (n!)
In SAS, this is computed as: EXP(-m) * (m**n) / FACT(n)

3. Using the recursive formula:
P(m,0) = exp(-m) and P(m,n) = P(m,(n-1)) * (m / n)


Using m = 5 and computing up to n = 30, I have found bigger differences between the methods than I would like to see. Methods 2 and 3 agree all the way to n = 30, but Method 1, using the SAS POISSON function, starts to diverge at n = 26.


CAN ANYONE PLEASE TELL ME WHICH METHOD IS THE MOST ACCURATE?



18 %let poimean = 5 ;
19
20
21 %let uplim = 30 ;

37 * method 2 ;
38 data _null_ ;
39 expm = exp(-&poimean.) ;
40 put expm= ;
41 do n = 0 to &uplim. ;
42 poi_n = expm * (&poimean.**n) / fact(n) ;
43 put n= poi_n= ;
44 end ;
45 run ;

expm=0.006737947
n=0 poi_n=0.006737947
n=1 poi_n=0.033689735
n=2 poi_n=0.0842243375
n=3 poi_n=0.1403738958
n=4 poi_n=0.1754673698
n=5 poi_n=0.1754673698
n=6 poi_n=0.1462228081
n=7 poi_n=0.104444863
n=8 poi_n=0.0652780393
n=9 poi_n=0.0362655774
n=10 poi_n=0.0181327887
n=11 poi_n=0.0082421767
n=12 poi_n=0.0034342403
n=13 poi_n=0.0013208616
n=14 poi_n=0.0004717363
n=15 poi_n=0.0001572454
n=16 poi_n=0.0000491392
n=17 poi_n=0.0000144527
n=18 poi_n=4.0146404E-6
n=19 poi_n=1.0564843E-6
n=20 poi_n=2.6412108E-7
n=21 poi_n=6.2885971E-8
n=22 poi_n=1.4292266E-8
n=23 poi_n=3.1070144E-9
n=24 poi_n=6.472947E-10
n=25 poi_n=1.294589E-10
n=26 poi_n=2.489595E-11
n=27 poi_n=4.610361E-12
n=28 poi_n=8.232787E-13
n=29 poi_n=1.419446E-13
n=30 poi_n=2.365743E-14
NOTE: DATA statement used (Total process time):
real time 0.15 seconds
cpu time 0.03 seconds


46
47 * method 3 ;
48 data _null_ ;
49 expm = exp(-&poimean.) ;
50 n = 0 ;
51 poi_n = expm ;
52 put n= poi_n= ;
53 do n = 1 to &uplim. ;
54 poi_n = poi_n * &poimean. / n ;
55 put n= poi_n= ;
56 end ;
57 run ;

n=0 poi_n=0.006737947
n=1 poi_n=0.033689735
n=2 poi_n=0.0842243375
n=3 poi_n=0.1403738958
n=4 poi_n=0.1754673698
n=5 poi_n=0.1754673698
n=6 poi_n=0.1462228081
n=7 poi_n=0.104444863
n=8 poi_n=0.0652780393
n=9 poi_n=0.0362655774
n=10 poi_n=0.0181327887
n=11 poi_n=0.0082421767
n=12 poi_n=0.0034342403
n=13 poi_n=0.0013208616
n=14 poi_n=0.0004717363
n=15 poi_n=0.0001572454
n=16 poi_n=0.0000491392
n=17 poi_n=0.0000144527
n=18 poi_n=4.0146404E-6
n=19 poi_n=1.0564843E-6
n=20 poi_n=2.6412108E-7
n=21 poi_n=6.2885971E-8
n=22 poi_n=1.4292266E-8
n=23 poi_n=3.1070144E-9
n=24 poi_n=6.472947E-10
n=25 poi_n=1.294589E-10
n=26 poi_n=2.489595E-11
n=27 poi_n=4.610361E-12
n=28 poi_n=8.232787E-13
n=29 poi_n=1.419446E-13
n=30 poi_n=2.365743E-14
NOTE: DATA statement used (Total process time):
real time 0.03 seconds
cpu time 0.03 seconds


58 * method 1 ;
59 data _null_ ;
60 n = 0 ;
61 poi_n = Poisson(&poimean.,n) ;
62 put n= poi_n= ;
63 nm1 = n ;
64 do n = 1 to &uplim. ;
65 poi_n = Poisson(&poimean.,n) - Poisson(&poimean.,nm1) ;
66 put n= poi_n= ;
67 nm1 = n ;
68 end ;
69 run ;

n=0 poi_n=0.006737947
n=1 poi_n=0.033689735
n=2 poi_n=0.0842243375
n=3 poi_n=0.1403738958
n=4 poi_n=0.1754673698
n=5 poi_n=0.1754673698
n=6 poi_n=0.1462228081
n=7 poi_n=0.104444863
n=8 poi_n=0.0652780393
n=9 poi_n=0.0362655774
n=10 poi_n=0.0181327887
n=11 poi_n=0.0082421767
n=12 poi_n=0.0034342403
n=13 poi_n=0.0013208616
n=14 poi_n=0.0004717363
n=15 poi_n=0.0001572454
n=16 poi_n=0.0000491392
n=17 poi_n=0.0000144527
n=18 poi_n=4.0146404E-6
n=19 poi_n=1.0564843E-6
n=20 poi_n=2.6412108E-7
n=21 poi_n=6.2885971E-8
n=22 poi_n=1.4292266E-8
n=23 poi_n=3.1070144E-9
n=24 poi_n=6.472947E-10
n=25 poi_n=1.294589E-10
n=26 poi_n=2.489597E-11
n=27 poi_n=4.610312E-12
n=28 poi_n=8.233414E-13
n=29 poi_n=1.418865E-13
n=30 poi_n=2.364775E-14


Thanks!

[Ksharp]

Hi.
I think the problem is n! .when n greater than 26 n! will be very very large to out the range of sas 's exact integer.Normally due to big n ,we will use Scott Formula
to approach n! .That is the reason why method 1 is different from 2 3.
In my opinion,Recommend to use method 1.
I am not quit sure my answer is right,it is just my opinion.


Ksharp