Re: Skellam distribution
- To: mathgroup at smc.vnet.net
- Subject: [mg102205] Re: [mg102152] Skellam distribution
- From: DrMajorBob <btreat1 at austin.rr.com>
- Date: Sat, 1 Aug 2009 04:00:33 -0400 (EDT)
- References: <200907310952.FAA19191@smc.vnet.net>
- Reply-to: drmajorbob at bigfoot.com
Here's the PDF from Mathematica,
f1 = PDF[d1 = PoissonDistribution[m1]];
f2 = PDF[d2 = PoissonDistribution[m2]];
pdf[k_] = Sum[f1[k + n] f2[n], {n, -k, Infinity}] // PowerExpand
E^(-m1 - m2) m1^(k/2) m2^(-k/2) BesselI[-k, 2 Sqrt[m1] Sqrt[m2]]
verified at
http://en.wikipedia.org/wiki/Skellam_distribution
and here's the CDF:
messy[k_] = Sum[pdf[n],{n,0,k}]
E^(-m1-m2) (-((Sqrt[m2] BesselI[1,2 Sqrt[m1]
Sqrt[m2]])/Sqrt[m1])+DifferenceRoot[Function[{\[FormalY],\[FormalN]},{m1
\[FormalY][\[FormalN]]+(-1-\[FormalN]-m1)
\[FormalY][1+\[FormalN]]+(1+\[FormalN]-m2) \[FormalY][2+\[FormalN]]+m2
\[FormalY][3+\[FormalN]]==0,\[FormalY][-1]==0,\[FormalY][0]==(Sqrt[m2]
BesselI[1,2 Sqrt[m1] Sqrt[m2]])/Sqrt[m1],\[FormalY][1]==BesselI[0,2
Sqrt[m1] Sqrt[m2]]+(Sqrt[m2] BesselI[1,2 Sqrt[m1]
Sqrt[m2]])/Sqrt[m1]}]][1+k])
That doesn't look nice at all, so I'd probably use
cdf[-1] = 0;
cdf[k_] := cdf[k] = cdf[k - 1] + pdf[k]
messy[10] - cdf[10] // FullSimplify
0
Bobby
On Fri, 31 Jul 2009 04:52:36 -0500, Peter Breitfeld <phbrf at t-online.de>
wrote:
>
> Is there an implementation of the Skellam distribution? If not, how
> would you implement a CDF for this distribution?
>
--
DrMajorBob at bigfoot.com