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!
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
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
Thanks @FreelanceReinh. This is much appreciated!! works great. Now onto the second half..the CDF.
Calling @Rick_SAS
Since FreelanceReinha
to calculate the cumulative probabilities, I added a line before the 'output' statement of the skellam datastep:
'cp+p;'
@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.
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!
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.