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