SAS function for the skellam distribution

[BillT]

I am looking for a function that will calculate the CDF of the skellam distribution. I didn't find the skellam distribution in the list of available SAS CDF functions. I realize I could code it but before I do that (and then test it...!), I thought I'd ask if the function is available.

Thanks!

[FreelanceReinh · Accepted solution]

Hello @BillT,

As far as I see, this distribution is not directly available in SAS. But, as you said, it's possible to code it.

For starters (and for curiosity), I've written this code for the PDF (and leave the cumulative summation for the CDF to you):

/* Compute pdf of the Skellam distribution */

%let mu1=5; /* parameter 1 */
%let mu2=2; /* parameter 2 */
%let c=7;   /* compute pdf for mean +/- c*sd */

data skellam(keep=x p);
m=&mu1-&mu2; /* mean */
t=&mu1+&mu2;
s=sqrt(t);   /* sd */
do x=floor(m-&c*s) to ceil(m+&c*s);
  p=exp(-t)*divide(&mu1,&mu2)**(x/2)*ibessel(abs(x),2*sqrt(&mu1*&mu2),0);
  output;
end;
run;

proc print data=skellam;
format p 16.14; /* e20.; */
run;

/* Simulate differences of independent Poisson random variables */

data sim;
call streaminit(27182818);
do i=1 to 1e6;
  x=rand('poisson',&mu1)-rand('poisson',&mu2);
  output;
end;
run;

proc means data=sim;
var x;
run; /* check mean and sd */

proc freq data=sim;
tables x / out=emp;
run;

/* Compare theoretical and empirical results */

data cmp(keep=x p f);
merge skellam
      emp(in=e);
by x;
if e then f=percent/100;
else f=0;
run;

proc sgplot data=cmp;
scatter x=x y=p;
scatter x=x y=f;
run;

The plot looks good. 🙂

Thanks for posting such an inspiring question!

Edit: I've used the formulas from Wikipedia: https://en.wikipedia.org/wiki/Skellam_distribution

[BillT]

Thanks @FreelanceReinh. This is much appreciated!! works great. Now onto the second half..the CDF.

[Ksharp]

Calling @Rick_SAS 

[Rick_SAS]

Since FreelanceReinhard has already provided the PDF, the only step left is to accumulate the PDF to form the CDF. As shown in the article "Create a customer PDF and CDF in SAS," there are two ways to proceed, depending upon how much precision you need.  In my article, I use the SAS/IML language, but the "easy" way can also be implemented in the DATA step by using the LAG function.

[BillT]

to calculate the cumulative probabilities, I added a line before the 'output' statement of the skellam datastep:
'cp+p;'

[FreelanceReinh]

---

@BillT wrote:

to calculate the cumulative probabilities, I added a line before the 'output' statement of the skellam datastep:
'cp+p;'

---

This should be fine. If in doubt, just increase the value of macro variable c. As soon as the largest computed cp values are equal to 1 (within the accuracy limits), you see that the omitted lower tail of the distribution (i.e., the sum of the p values for x<floor(m-&c*s)) was negligible.