Here's the notebook for Poisson Distributions. Just select it then Copy and Paste to your Mathematica notebook. For convenient presentation it's a single notebook, but for easier use divide it into cells. If you do so, be sure to also load packages (<<) needed. To test your output compare with the figures in this web.
  The last few lines of this notebook use the function Display[ ] to save each of the graphics as GIF files.

 

(*   Poisson Probability Distribution

     To accompany "Poisson Distributions"

     William J. Thompson,

     Computing in Science and Engineering,

     May/June 2001 *)

 

Off[TagSetDelayed::write,General::spell1,Display::Pserr]

<<Statistics`DiscreteDistributions`

<<Graphics`Graphics3D`

<<Graphics`Graphics`

 

(* Poisson Distribution table *)

PoissDistTable[m_,kmin_,kmax_] := Table[

  PDF[PoissonDistribution[m],k],{k,kmin,kmax}];

 

(* Poisson Distribution: 3D views *)

(* Figure 1(a) *)

Poisson3D =

  BarChart3D[

  {PoissDistTable[4.1,1,26],PoissDistTable[9.1,1,26],

  PoissDistTable[16.1,1,26]},Frame->False,

  Ticks->{None,{5,10,15,20,25},None},

  DefaultFont->{"Helvetica",14},

  AxesEdge->{{1,-1},{1,-1},Automatic},

  AxesLabel->{"","  k"," "},ViewPoint->{2.5,-0.5,1}];

    

(* Poisson Distribution; 2D views *)

PoissonPlot[k_] := Plot[PDF[PoissonDistribution[m],k],

  {m,0,20},PlotRange->All,Ticks->{Automatic,None},

  DefaultFont->{"Helvetica",14},AxesLabel->{"m ",""},

  DisplayFunction->Identity];

    

(* Figure 1(b) *)

ShowPoisson = Show[PoissonPlot[2],PoissonPlot[4],

  PoissonPlot[9],DisplayFunction->$DisplayFunction];

 

(* Comparison of distributions *)

(* Binomial Distribution table *)

BinomDistTable[N_,p_] := Table[

  PDF[BinomialDistribution[N,p],k],{k,1,20}];

BinomialPoisson[N_,m_] := BarChart[

  PoissDistTable[m,1,20],BinomDistTable[N,m/N],

  PlotRange->All,Ticks->{Automatic,None},

  AxesLabel->{"k",""},DefaultFont->{"Helvetica",30},

  PlotLabel->"p = "<>ToString[m/N]<>"",

  BarSpacing->0.8,BarGroupSpacing->0.3,BarEdges->False,

  BarStyle->{Hue[0.99,0.9,0.9],Hue[0.5,0.6,0.95]}];

 

(* Figure 2 *)

BP1041 = BinomialPoisson[10,4.1];

BP1091 = BinomialPoisson[10,9.1];

BP10041 = BinomialPoisson[100,4.1];

BP10091 = BinomialPoisson[100,9.1];

 

(* Gauss distribution *)

GaussProb[m_,x_] := (1/(Sqrt[2*N[Pi]*m]))*

  Exp[-(x-m)^2/(2*m)];

GaussPlot[m_]:=

  Plot[GaussProb[m,x],{x,1,20},PlotPoints->100,

  Ticks->{Automatic,None},

  AxesLabel->{"k",""},DefaultFont->{"Helvetica",20},

  PlotStyle->{Thickness[0.008]},DisplayFunction->Identity];

PoissonPlot[m_] := BarChart[PoissDistTable[m,1,20],

  PlotRange->All,Ticks->{Automatic,None},

  AxesLabel->{"k",""},DefaultFont->{"Helvetica",20},

  PlotLabel->"m = "<>ToString[m]<>"",

  BarStyle->{Hue[0.99,0.9,0.9]},BarEdges->False,

  DisplayFunction->Identity]

GaussPoisson[m_]:=

  Show[PoissonPlot[m],GaussPlot[m],

  DisplayFunction->$DisplayFunction];

   

(* Figure 3(a) *)

GP41 = GaussPoisson[4.1];

(* Figure 3(b) *)

GP91 = GaussPoisson[9.1];

 

(* RMS fractional error in Gauss distribution

   compared to Poisson distribution *)

(* Gaussian averaged over k-1/2 to k+1/2 *)

PGAver[m_,k_] :=

  (Erf[(k-m+1/2)/Sqrt[2m]]-

   Erf[(k-m-1/2)/Sqrt[2m]])/2

RMS[m_] :=

  Sqrt[NSum[(PGAver[m,k]

  -PDF[PoissonDistribution[m],k])^2,{k,0,20,1}]];

 

(* Figure 4 *)

RMSPlot =

  Plot[RMS[m],{m,0.5,20},PlotPoints->200,

  Ticks->{{5,10,15,20},Automatic},

  AxesLabel->{" m","RMS (m)"},DefaultFont->{"Helvetica",14},

  PlotStyle->{Thickness[0.004]}];

 

(* Distribution of subtracted Poisson variables *)

(* Distribution function *)

PoissSub[m1_,m2_,d_] := If[

  d>=0, Exp[-(m1+m2)]*

  (m1/m2)^(d/2)*BesselI[d,2*Sqrt[m1*m2]],

   PoissSub[m2,m1,-d]]

(* Table of subtracted distribution *)

PoissSubTable[m1_,m2_] := Table[

  PoissSub[m1,m2,d],{d,-10,10}]

 

(* Define bar chart *)

PoissSubBar[m1_,m2_] := BarChart[

  PoissSubTable[m1,m2],PlotRange->All,

  BarStyle->{Hue[0.35,0.9,1]},

  BarEdges->False ,Axes->{True,False},

  AxesLabel->{"d",None},Ticks->False,

  DefaultFont->{"Helvetica",20}]

 

(* Figure 5 *)

PoissSubPic11 = PoissSubBar[1.1,1.1];

PoissSubPic14 = PoissSubBar[1.1,4.1];

PoissSubPic41 = PoissSubBar[4.1,1.1];

PoissSubPic44 = PoissSubBar[4.1,4.1];

 

(* Moments *)

(* Defined about the origin *)

PoissMomZero[m1_,m2_,n_] :

  NSum[d^n*PoissSub[m1,m2,d],{d,-10,10}]

(* Defined about the mean *)

PoissMomMean[m1_,m2_,n_] :=

  NSum[(d-(m1-m2))^n*PoissSub[m1,m2,d],

  {d,-10,10}]

(* Example in Figure 5 *)

m1 = 1.1; m2 = 4.1;

SubMean = m1-m2

PoissMomZero[m1,m2,1]

PoissMomMean[m1,m2,1]

SubVariance = m1+m2

PoissMomZero[m1,m2,2]

PoissMomMean[m1,m2,2]

ThirdMoment = m1-m2

PoissMomZero[m1,m2,3]

PoissMomMean[m1,m2,3]

 

(* Poisson-distributed random numbers *)

(* Algorithm for a random number from

   a Poisson distribution with mean = m *)

PoissRandom[m_] := (

  RandomProd = 1;

  Do[RandomProd *= Random[];

    iKeep = i;

    If[RandomProd<Exp[-m],Break[]],

    {i,1,Max[10,2*m]}];

    Return[iKeep-1])

(* Define table of Nran Poisson random numbers *)

PoissRanList[m_,Nran_] := Table[

  PoissRandom[m],{i,1,Nran}]

 

m = 4.1;

Nran = 100;

 

(* Table of how many of each k, up to 3m - 1 *)

PoissRandomDist =

Table[Count[Evaluate[PoissRanList[m,Nran]],k-1],

  {k,3*m}];

(* Total for normalization *)

NormRan = Apply[Plus,PoissRandomDist];

PoissFormulaTable = Table[

  NormRan*Exp[-m]*m^k/k!,{k,0,3*m-1}];

(* Bar chart *)

(* Figure 6(a). For Figure 6(b) set Nran = 1000 above.

   Note: PoissRandomDist changes from run to run *)

RanForm = BarChart[

  PoissRandomDist,PoissFormulaTable,

  PlotRange->All,Axes->{True,False},Ticks->None,

  AxesLabel->{"k",""},BarEdges->False,

  DefaultFont->{"Helvetica",20},

  BarSpacing->0.8,BarGroupSpacing->0.3,

  BarStyle->{Hue[0.7,0.8,0.8],Hue[0.99,0.9,0.9]}];

 

 

(* To save graphics as GIFs.

      The GIFs will likely be in your Mathematica folder *)

 

Off[Display::pserr];

(* Figure 1(a) *)

Display["Poisson3D.GIF", Poisson3D, "GIF"];

(* Figure 1(b) *)

Display["ShowPoisson.GIF", ShowPoisson, "GIF"];

(* Figure 2 *)

Display["BP1041.GIF", BP1041, "GIF"];

Display["BP1091.GIF", BP1091, "GIF"];

Display["BP10041.GIF", BP10041, "GIF"];

Display["BP10091.GIF", BP10091, "GIF"];

(* Figure 3 *)

Display["GP41.GIF", GP41, "GIF"];

Display["GP91.GIF", GP91, "GIF"];

(* Figure 4 *)

Display["RMSPlot.GIF", RMSPlot, "GIF"];

(* Figure 5 *)

Display["PoissSubPic11.GIF", PoissSubPic11, "GIF"];

Display["PoissSubPic14.GIF", PoissSubPic14, "GIF"];

Display["PoissSubPic41.GIF", PoissSubPic41, "GIF"];

Display["PoissSubPic44.GIF", PoissSubPic44, "GIF"];

(* Figure 6(a) *)

Display["RanForm.GIF", RanForm, "GIF"];