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

Re: Skellam distribution

• To: mathgroup at smc.vnet.net
• Subject: [mg102222] Re: [mg102152] Skellam distribution
• From: Bob Hanlon <hanlonr at cox.net>
• Date: Sat, 1 Aug 2009 04:04:02 -0400 (EDT)
• Reply-to: hanlonr at cox.net

Clear[SkellamDistribution];

SkellamDistribution::usage =
  "SkellamDistribution[a1,a2] is the Skellam discrete probability \
distribution with parameters a1 and a2. \
http://en.wikipedia.org/wiki/Skellam_distribution";;

SkellamDistribution /: PDF[
   SkellamDistribution[a1_, a2_], k_] :=
  Exp[-(a1 + a2)]*(a1/a2)^(k/2)*
   BesselI[k, 2 Sqrt[a1*a2]];

SkellamDistribution /: CDF[
   SkellamDistribution[a1_, a2_], k_] :=
  Sum[PDF[
    SkellamDistribution[a1, a2], x],
   {x, -Infinity, k}];

Off[SequenceLimit::seqlim];

{#, N[CDF[SkellamDistribution[4, 2], #]]} & /@ {-10, 0, 5, 10, 15, 20, 21, =
22,
   23, 24}

{{-10, 1.7071667532058565*^-6}, {0, 0.2700394539486425},
   {5, 0.9224208128081963}, {10, 0.9992781041223845},
   {15, 0.9999989537444384}, {20, 0.9995347398065093},
   {21, 0.999692263837724}, {22, 1.0025003178112426},
   {23, 0.00034632170838100936},
   {24, 0.00025879506570215534}}

Although this is a precision problem; just increasing the precision doesn't=
 help much.

Off[N::meprec];

Block[{$MaxExtraPrecision = 200},
 {#, N[CDF[SkellamDistribution[4, 2], #], 50]} & /@
  {20, 21, 22, 23, 24}]

{{20,
     0.99953473980650827322371799398032236022`2.744357869581\
   4312},
   {21,
     0.99969226383772326989075589292344256806`3.628450699957\
   981},
   {22,
     1.002500317811242178642646033555563141`2.41501144500697\
   86},
   {23,
     0.00034632170838105499545748845259004187`0.890238973328\
   0173},
   {24,
     0.00025879506570205721307174588171233687`1.237949261004\
   8972}}

Instead, aggregating the sum differently for values above the mean appears =
to work well.

SkellamDistribution /: CDF[
   SkellamDistribution[a1_, a2_], k_Integer] :=
 
  PDF[SkellamDistribution[a1, a2], a1 - a2] +
    Sum[
     PDF[SkellamDistribution[a1, a2], a1 - a2 + n] +
      PDF[SkellamDistribution[a1, a2], a1 - a2 - n],
     {n, 1, k - (a1 - a2)}] +
    Sum[
     PDF[SkellamDistribution[a1, a2], n],
     {n, -Infinity, 2 (a1 - a2) - k - 1}] /; k > (a1 - a2);

ListPlot[{#, CDF[SkellamDistribution[4, 2], #]} & /@
  Range[-10, 40],
 Filling -> Axis,
 Frame -> True,
 Axes -> False]

An alternative grouping could be investigated if this not work for a partic=
ular problem.


Bob Hanlon

---- Peter Breitfeld <phbrf at t-online.de> wrote:

=============
Thank you very much Bob!

All of these work fine, except the CDF. I did implement it the same 
way you did, but not using UpValues. But using the CDF either way look 
at this output, which gives big numerical errors

In=({#, CDF[SkellamDistribution[4, 2], #]} & /@
     {-10, 0, 5, 10, 15, 20, 21, 22, 23, 24}
    ) // N

SequenceLimit::seqlim :
The general form of the sequence could not be
determined, and the result may be incorrect.

Out={{-10., 1.70717*10^-6}, {0., 0.270039}, {5., 0.922421},
      {10., 0.999278}, {15., 0.999999}, {20., 0.999535}, {21., 
0.999692},
      {22., 1.0025}, {23., 0.000346322}, {24., 0.000258795}}

It seems kind of chaotic behaviour takes place for bigger values. In 
this example all seems ok up to 21, but 22 gives a value greater than 
1, which can't be and the folliwing values switch down to about Zero.

I tried to give the Sum used in CDF options like Method or 
Regularization, but they give errors, because Mathematica replaces 
N[Sum[...]] with NSum, which doesn't know these options.

Am 31. Jul 2009 um 14:06 schrieb Bob Hanlon:

>
> Clear[SkellamDistribution];
>
> SkellamDistribution::usage =
>  "SkellamDistribution[a1,a2] is the Skellam discrete probability \
> distribution with parameters a1 and a2. \
> http://en.wikipedia.org/wiki/Skellam_distribution";;
>
> SkellamDistribution /: PDF[
>   SkellamDistribution[a1_, a2_], k_] :=
>
>  Exp[-(a1 + a2)]*(a1/a2)^(k/2)*
>   BesselI[k, 2 Sqrt[a1*a2]];
>
> SkellamDistribution /: mgf[
>   SkellamDistribution[a1_, a2_], t_] :=
>
>  Exp[-(a1 + a2) + a1*Exp[t] + a2*Exp[-t]];
>
> SkellamDistribution /: CharacteristicFunction[
>   SkellamDistribution[a1_, a2_], t_] :=
>
>  Exp[-(a1 + a2) + a1*Exp[I*t] + a2*Exp[-I*t]];
>
> SkellamDistribution /: Moment[
>   SkellamDistribution[a1_, a2_], n_Integer] :=
>  Module[{t},
>   D[mgf[SkellamDistribution[a1, a2], t], {t, n}] /.
>    t -> 0];
>
> SkellamDistribution /: CDF[
>   SkellamDistribution[a1_, a2_], k_] :=
>  Sum[PDF[
>    SkellamDistribution[a1, a2], x],
>   {x, -Infinity, k}];
>
> SkellamDistribution /: Mean[
>   SkellamDistribution[a1_, a2_]] :=
>  (a1 - a2);
>
> SkellamDistribution /: Variance[
>   SkellamDistribution[a1_, a2_]] :=
>  (a1 + a2);
>
> SkellamDistribution /: StandardDeviation[
>   SkellamDistribution[a1_, a2_]] :=
>  Sqrt[a1 + a2];
>
> Table[{n,
>  Moment[SkellamDistribution[a1, a2], n]},
> {n, 0, 2}]
>
> {{0, 1}, {1, a1 - a2}, {2, (a1 - a2)^2 + a1 + a2}}
>
>
> Bob Hanlon
>
> ---- 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?
>
> --
> _________________________________________________________________
> Peter Breitfeld, Bad Saulgau, Germany -- http://www.pBreitfeld.de
>
>

Gru=C3=9F Peter
--
__________________________________________________________
Peter Breitfeld                 | http://www.pBreitfeld.de
Kreuzgasse 4, 88348 Bad Saulgau | PGP/GPG Key gibt's hier