topic Re: SAS function for the skellam distribution in Statistical Procedures

SAS function for the skellam distribution

BillT — Sat, 11 Jan 2020 16:53:03 GMT

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!

Re: SAS function for the skellam distribution

FreelanceReinh — Sat, 11 Jan 2020 19:51:17 GMT

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

Re: SAS function for the skellam distribution

Ksharp — Sun, 12 Jan 2020 11:03:44 GMT

Calling @Rick_SAS

Re: SAS function for the skellam distribution

Rick_SAS — Sun, 12 Jan 2020 11:28:51 GMT

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.

Re: SAS function for the skellam distribution

BillT — Tue, 14 Jan 2020 00:16:14 GMT

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

Re: SAS function for the skellam distribution

BillT — Wed, 15 Jan 2020 16:22:42 GMT

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

Re: SAS function for the skellam distribution

FreelanceReinh — Wed, 15 Jan 2020 17:12:23 GMT

@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.