/* Calculate theoretical Negative Binomial PMF based on fitted k and p */;
***In SAS, m is number of failures, p is probability of success, k is number of failures;
***k = 1/dispersion; ****p = 1/(1+exp(intercept)*dispersion);
data TheoreticalNegBin1;
   k=0.1;
   p=0.5;
   mincount=0;
   maxcount=20;
   do count=minCount to maxCount;
      NegBin1=pdf('NegBinomial', count, p, k);
      output;
   end;
run; 

data TheoreticalNegBin2;
   k=1;
   p=0.5;
   mincount=0;
   maxcount=20;
   do count=minCount to maxCount;
      NegBin2=pdf('NegBinomial', count, p, k);
      output;
   end;
run; 

data TheoreticalNegBin3;
   k=2;
   p=0.5;
   mincount=0;
   maxcount=20;
   do count=minCount to maxCount;
      NegBin3=pdf('NegBinomial', count, p, k);
      output;
   end;
run; 

data TheoreticalNegBin4;
   k=3;
   p=0.5;
   mincount=0;
   maxcount=20;
   do count=minCount to maxCount;
      NegBin4=pdf('NegBinomial', count, p, k);
      output;
   end;
run;

data TheoreticalNegBin5;
   k=10;
   p=0.5;
   mincount=0;
   maxcount=20;
   do count=minCount to maxCount;
      NegBin5=pdf('NegBinomial', count, p, k);
      output;
   end;
run;

data TheoreticalNegbinAll;
 merge TheoreticalNegBin1 TheoreticalNegBin2 TheoreticalNegBin3 TheoreticalNegBin4 TheoreticalNegBin5;
 by count;
 run;
 
title 'Theoretical distributions';
proc sgplot data=TheoreticalNegbinAll noautolegend;
   series x=count y=Negbin1 / markers markerattrs=(symbol=circlefilled color=red) 
                         lineattrs=(color=red)legendlabel='Negative Binomial k=0.1';
   series x=count y=Negbin2 / markers markerattrs=(symbol=circlefilled color=yellow) 
                         lineattrs=(color=yellow)legendlabel='Negative Binomial k=1'; 
   series x=count y=Negbin3 / markers markerattrs=(symbol=circlefilled color=green) 
                         lineattrs=(color=green)legendlabel='Negative Binomial k=2';
   series x=count y=Negbin4 / markers markerattrs=(symbol=circlefilled color=blue) 
                         lineattrs=(color=blue)legendlabel='Negative Binomial k=3'; 
   series x=count y=Negbin5 / markers markerattrs=(symbol=circlefilled color=violet) 
                         lineattrs=(color=violet)legendlabel='Negative Binomial k=10'; 
   xaxis label = 'COUNT' fitpolicy=THIN;
   yaxis label = 'PROBABILTY MASS FUNCTION';
   keylegend / location=inside position=NE across=1;
run;
title;