BookmarkSubscribeRSS Feed
🔒 This topic is solved and locked. Need further help from the community? Please sign in and ask a new question.
BillT
Fluorite | Level 6

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!

1 ACCEPTED SOLUTION

Accepted Solutions
FreelanceReinh
Jade | Level 19

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

View solution in original post

6 REPLIES 6
FreelanceReinh
Jade | Level 19

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
Fluorite | Level 6

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

Rick_SAS
SAS Super FREQ

Since 

BillT
Fluorite | Level 6

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

 

 

FreelanceReinh
Jade | Level 19

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

SAS Innovate 2025: Register Now

Registration is now open for SAS Innovate 2025 , our biggest and most exciting global event of the year! Join us in Orlando, FL, May 6-9.
Sign up by Dec. 31 to get the 2024 rate of just $495.
Register now!

What is ANOVA?

ANOVA, or Analysis Of Variance, is used to compare the averages or means of two or more populations to better understand how they differ. Watch this tutorial for more.

Find more tutorials on the SAS Users YouTube channel.

Discussion stats
  • 6 replies
  • 1230 views
  • 4 likes
  • 4 in conversation