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)

```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

```

• Prev by Date: Re: Finding the Position of Elements in a List that
• Next by Date: Re: Mapping GPS data
• Previous by thread: Re: Skellam distribution
• Next by thread: Re: How to define a list of functions/variables