MathGroup Archive 2009

[Date Index] [Thread Index] [Author Index]

Search the Archive

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


  • Prev by Date: Re: Skellam distribution
  • Next by Date: Re: Finding the Position of Elements in a List that Contain
  • Previous by thread: Re: Skellam distribution
  • Next by thread: Re: Skellam distribution