Re: Skellam distribution (solved)

• To: mathgroup at smc.vnet.net
• Subject: [mg102215] Re: Skellam distribution (solved)
• From: Peter Breitfeld <phbrf at t-online.de>
• Date: Sat, 1 Aug 2009 04:02:37 -0400 (EDT)
• References: <h4uesu\$ink\$1@smc.vnet.net>

```Peter Breitfeld wrote:

> Is there an implementation of the Skellam distribution? If not, how
> would you implement a CDF for this distribution?
>
> --
> _________________________________________________________________
> Peter Breitfeld, Bad Saulgau, Germany -- http://www.pBreitfeld.de
>

With the help of Bob Hanlon, Darren Glosemeyer and DrMajor Bob (thank
you very much!) I could solve the main problem I had with this: getting
the numerical errors away.

My first solution was to write:

skellamCDF[l_,m_,k_]:=
NSum[E^(-l-m) (l/m)^(k/2) BesselI[x, 2 Sqrt[l m]],{x,-Infinity,k}];

But this gave heavy rounding errors e.g.:

In= skellanCDF[2.5, 3, #] & /@ {-10, 0, 5, 10, 15, 20, 21, 22, 23, 24}

Many messages <SequenceLimit::seqlim>

Out= {0.000172756, 0.667993, 0.994528, 0.999994, 1., 1.59882,
0.00017627, -0.0000221875, 0.0001553, 0.000321263}

As you can see, these are unusable results.

My implementation now  seems to work fine, and, thanks to Bob,
behaves like other build in distributions:

SkellamDistribution /: PDF[
SkellamDistribution[l_Positive, m_Positive], k_Integer] :=
E^-(l + m) (l/m)^(k/2) BesselI[k, 2 Sqrt[l m]];

SkellamDistribution /: CharacteristicFunction[
SkellamDistribution[l_?Positive, m_?Positive], t_?NumericQ] :=
Exp[-(l + m) + l Exp[I t] + m Exp[-I t]];

SkellamDistribution /: CDF[
SkellamDistribution[l_?Positive, m_?Positive], k_?NumericQ] :=
If[k > l - m, (* l-m is the Mean *)
1 - NSum[PDF[SkellamDistribution[l, m], x],
{x, Ceiling[k + 1], Infinity}],
NSum[PDF[SkellamDistribution[l, m], x],
{x, -Infinity, Floor[k]}]];

SkellamDistribution /: Mean[
SkellamDistribution[l_?Positive, m_?Positive]] := (l - m);

SkellamDistribution /: Variance[
SkellamDistribution[l_?Positive, m_?Positive]] := (l + m);

SkellamDistribution /: StandardDeviation[
SkellamDistribution[l_?Positive, m_?Positive]] := Sqrt[l + m];

With this CDF I get:

In= {#, CDF[SkellamDistribution[2.5, 3], #]} & /@
{-10, 0, 5, 10, 15, 20, 21, 22, 23, 24}

Out= {{-10, 0.000172756}, {0, 0.667993}, {5, 0.994528}, {10, 0.999994},
{15, 1.}, {20, 1.}, {21, 1.}, {22, 1.}, {23, 1.}, {24, 1.}}

which seems to be very reasonable.

--
_________________________________________________________________
Peter Breitfeld, Bad Saulgau, Germany -- http://www.pBreitfeld.de

```

• Prev by Date: Re: Finding the Position of Elements in a List that
• Next by Date: Re: Problem in plotting Bifurcation Diagram (ListPlot with
• Previous by thread: Re: Mapping GPS data
• Next by thread: Re: Problem in plotting Bifurcation Diagram (ListPlot with