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

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