Discrete Lerch distribution

This package extends Mathematica standard statistical functions to include the discrete Lerch distribution and itrs special cases, the Zipf, Zipf-Mandelbrot, and Good distributions.  In addition, Lerch distribution and its special cases can be fitted to data with the method of moments, maximum likelihood, and by minimzation of Pearson's Chi-squared statistic.

Author:  Sergej V. Aksenov

Department of Microbiology and Immunology
University of Michigan
Ann Arbor, MI  48109, USA
http://aksenov.freeshell.org/

Last updated 17 June 2002

Installation and initialization

The package file LerchDistribution.m should be located inside the LerchDistribution directory.  
The LerchDistribution directory should be put in one of the canonical places for AddOns.   
Two recommended places are first inside the Mathematica distribution in:

ToFileName[ {$TopDirectory,  "AddOns",  "Applications"}]

a second in the preferences directory:

\(ToFileName[\ {$PreferencesDirectory, \ "AddOns", \ "Applications"}]\)

When one of these locations is used there is no need to set $Path and the loading instructions are identical for all versions of Mathematica.

The  LerchDistribution` package can  be loaded with

<< LerchDistribution`LerchDistribution`

Load the Lerch Phi by CNCT package (get it at http://aksenov.freeshell.org/)

<< lerchphi`lerchphi`

Introduction

Lerch distribution with three parameters z, s, and v are defined by the following probability mass function

(Subscript[p, x]) = (1/c) ((z^x)/(((v + x))^s)),  x = 0, 1, 
    2,  ...

where

c = Φ (z, s, v)

is the series function of the Lerch distribution (Lerch's transcendent) that insures that probabilities always sum up to unity.

For properties, functionals and estimation methods for the Lerch distribution see article
"Parameter estimation and random number generation from a Zipf-related Lerch distribution", S.V. Aksenov, M.A. Savageau,  submitted to Computational Statistics and Data Analysis, 2001.  Preprint MPS: Statistics/0206008 is available at http://www.mathpreprints.com/Statistics/0206008/

Special cases of the Lerch distribution, the Zipf (v=0, z=1), Zipf-Mandelbrot (z=1), and Good (v=0) distributions, and their properties, are discussed in the article
"Application of the Combined Nonlinear--Condensation Transformation to Problems in Statistical Analysis and Theoretical Physics", S.V. Aksenov, M.A. Savageau, U.D. Jentschura, J. Becher, G. Soff, P.J. Mohr, submitted to Computer Physics Communications, 2002.

Examples

<< Graphics`Graphics`

<< Statistics`ContinuousDistributions`

Auxilliary macros

pairs[s_List, n_Integer] := 
        Module[{i, res}, 
          res = {}; 
Do[
            AppendTo[res, Take[s, {i, i + n - 1}]], {i, 1, 
              Length[s] - n + 1}]; 
Return[
            res]];

    
putfreq[ff_List, cc_List, cuts_List] := 
        Module[{prs, cat1, inp, inp2, pos, fs}, 
          prs = pairs[cuts, 2]; 
cat1 = 
            Map[(inp = #; 
                  Select[cc, (# ≥ inp[[1]] && # < 
                            inp[[2]]) &]) &, 
              Drop[prs, -1]]; 
AppendTo[cat1, 
            Flatten[Map[(inp = #; 
Select[
                      cc, (# ≥ inp[[1]] && # ≤ 
                              inp[[2]]) &]) &, {Last[
                    prs]}]]]; 
pos = 
            Map[(inp2 = #; Flatten[Map[Position[cc, #] &, inp2]]) &, 
              cat1]; 
fs = Map[ff[[#]] &, pos]; 
          Map[Apply[Plus, #] &, fs]];

Example 1 (PDF, CDF)

Example 1 demonstrates the use of standard statistical functions with the Lerch distribution.

?LerchDistribution

Domain[LerchDistribution[0.7, 0.5, 0.5]]

Domain[LerchDistribution[0.7, 0.5, 0.5, Interval[{1, 100}]]]

DomainQ[LerchDistribution[0.7, 0.5, 0.5], 2]

DomainQ[LerchDistribution[0.7, 0.5, 0.5], -10]

DomainQ[LerchDistribution[0.7, 0.5, 0.5, Interval[{1, 100}]], 
      200]

PDF[LerchDistribution[0.7, 0.5, 0.5], 0]

PDF[LerchDistribution[0.7, 0.5, 0.5], Infinity]

PDF[LerchDistribution[0.7, 0.5, 0.5], -10]

Sum[
      PDF[LerchDistribution[0.99, 0.5, 0.5], x], {x, 0, Infinity}]

CDF[LerchDistribution[0.7, 0.5, 0.5], 0]

CDF[LerchDistribution[0.7, 0.5, 0.5], Infinity]

CDF[LerchDistribution[0.7, 0.5, 0.5], -10]

Example 1a (PDF, CDF  for Zipf etc.)

Example 1a demonstrates the use of standard statistical functions with special cases of the Lerch distribution.

Zipf (zeta)

Domain[
      LerchDistribution[1, 1.5, 0, Interval[{1, Infinity}]]]

DomainQ[LerchDistribution[1, 1.5, 0, Interval[{1, Infinity}]], 
      2]

DomainQ[LerchDistribution[1, 1.5, 0, Interval[{1, Infinity}]], 
      0]

PDF[LerchDistribution[1, 1.5, 0, Interval[{1, Infinity}]], 
      1]

Sum[
      PDF[LerchDistribution[1, 1.5, 0, Interval[{1, Infinity}]], x], {x, 1, 
        Infinity}]

CDF[LerchDistribution[1, 1.5, 0, Interval[{1, Infinity}]], 
      1]

CDF[LerchDistribution[1, 1.5, 0, Interval[{1, Infinity}]], 
      Infinity]

Zipf-Mandelbrot (generalized zeta)

Domain[LerchDistribution[1, 1.5, 2]]

DomainQ[LerchDistribution[1, 1.5, 2], 2]

DomainQ[LerchDistribution[1, 0.5, 2], -1]

PDF[LerchDistribution[1, 1.5, 2], 0]

Sum[PDF[LerchDistribution[1, 1.5, 2], x], {x, 0, Infinity}]

CDF[LerchDistribution[1, 1.5, 2], 0]

CDF[LerchDistribution[1, 1.5, 2], Infinity]

CDF[LerchDistribution[1, 1.5, 2], -1]

Good (Jonquiere)

Domain[
      LerchDistribution[0.99, 0.5, 0, Interval[{1, Infinity}]]]

DomainQ[LerchDistribution[0.99, 0.5, 0, Interval[{1, Infinity}]], 
      2]

DomainQ[LerchDistribution[0.99, 0.5, 0, Interval[{1, Infinity}]], 
      0]

PDF[LerchDistribution[0.99, 0.5, 0, Interval[{1, Infinity}]], 
      1]

Sum[
      PDF[LerchDistribution[0.99, 0.5, 0, Interval[{1, Infinity}]], x], {x, 
        1, Infinity}]

CDF[LerchDistribution[0.99, 0.5, 0, Interval[{1, Infinity}]], 
      1]

CDF[LerchDistribution[0.99, 0.5, 0, Interval[{1, Infinity}]], 
      Infinity]

CDF[LerchDistribution[0.99, 0.5, 0, Interval[{1, Infinity}]], 
      0]

Example 2 (Moments)

Example 2 demonstrates the use of standard functions of descriptive statistics with the Lerch distribution.

?UncorrectedMoment

UncorrectedMoment[LerchDistribution[0.7, 0.5, 0.5], 1]

Mean[LerchDistribution[0.7, 0.5, 0.5]]

ExpectedValue[# &, LerchDistribution[0.7, 0.5, 0.5]]

N[ExpectedValue[x, LerchDistribution[0.7, 0.5, 0.5], x]]

Variance[LerchDistribution[0.7, 0.5, 0.5]]

StandardDeviation[LerchDistribution[0.7, 0.5, 0.5]]

?CentralMoment

CentralMoment[LerchDistribution[0.7, 0.5, 0.5], 2]

ExpectedValue[(# - 1.2995094289063234`)^2 &, 
      LerchDistribution[0.7, 0.5, 0.5]]

?DescendingFactorialMoment

DescendingFactorialMoment[LerchDistribution[0.7, 0.5, 0.5], 
      4]

Skewness[LerchDistribution[0.7, 0.5, 0.5]]

Kurtosis[LerchDistribution[0.7, 0.5, 0.5]]

CharacteristicFunction[LerchDistribution[0.7, 0.5, 0.5], t]

Quantile[LerchDistribution[0.7, 0.5, 0.5], 0.9999]

Check with CDF

CDF[LerchDistribution[0.7, 0.5, 0.5], 20]
CDF[LerchDistribution[0.7, 0.5, 0.5], 21]

Quantile[LerchDistribution[0.7, 0.5, 0.5], 0.51]

Check with CDF

CDF[LerchDistribution[0.7, 0.5, 0.5], 0]
CDF[LerchDistribution[0.7, 0.5, 0.5], 1]

Example 2a (Moments for Zipf etc.)

Example 2a demonstrates the use of standard functions of descriptive statistics with special cases of the Lerch distribution.

Zipf (zeta)

UncorrectedMoment[
      LerchDistribution[1, 2.5, 0, Interval[{1, Infinity}]], 1]

Mean[LerchDistribution[1, 2.5, 0, Interval[{1, Infinity}]]]

Variance[
      LerchDistribution[1, 3.5, 0, Interval[{1, Infinity}]]]

StandardDeviation[
      LerchDistribution[1, 3.5, 0, Interval[{1, Infinity}]]]

CentralMoment[LerchDistribution[1, 3.5, 0, Interval[{1, Infinity}]], 
      2]

Zipf-Mandelbrot (generalized zeta)

UncorrectedMoment[LerchDistribution[1, 2.9, 1], 1]

Mean[LerchDistribution[1, 2.9, 1]]

Variance[LerchDistribution[1, 3.9, 1]]

StandardDeviation[LerchDistribution[1, 3.9, 1]]

CentralMoment[LerchDistribution[1, 3.9, 1], 2]

Good (Jonquiere)

UncorrectedMoment[
      LerchDistribution[0.9, 2.5, 0, Interval[{1, Infinity}]], 1]

Mean[
      LerchDistribution[0.9, 2.5, 0, Interval[{1, Infinity}]]]

Variance[
      LerchDistribution[0.9, 3.5, 0, Interval[{1, Infinity}]]]

StandardDeviation[
      LerchDistribution[0.9, 3.5, 0, Interval[{1, Infinity}]]]

CentralMoment[LerchDistribution[0.9, 3.5, 0, Interval[{1, Infinity}]], 
      2]

Example 3 (Random number generation)

Example 3 demonstrates random number generation from the Lerch distribution.

Random[LerchDistribution[0.7, 0.5, 0.5]]

SeedRandom[1121]

Timing[data = 
          RandomArray[LerchDistribution[0.7, 0.5, 0.5], 1000];]

ListPlot[Map[Reverse, Frequencies[data]], 
        PlotJoined -> True];

freq1 = Transpose[Frequencies[data]][[1]]

classes1 = Transpose[Frequencies[data]][[2]]

population = 
      Map[PDF[LerchDistribution[0.7, 0.5, 0.5], #] &, classes1]

datahist = 
        Histogram[data, DisplayFunction -> Identity];

popul = 
        ListPlot[Transpose[{classes1,  Length[data] population}], 
          PlotJoined -> True, PlotStyle -> {Hue[0.65]}, 
          DisplayFunction -> Identity];

Show[datahist, popul, 
        DisplayFunction -> $DisplayFunction];

Example 4 (MM and MLE)

Example 4 demonstrates how to fit Lerch distribution to data using method of moments and maximum likelihood.

NOTE: please generate data in Example 3  first!!

Initial guesses

Fit Zeta distribution to ratio of second frequency (class x=1)

sest = Log[2, freq1[[2]]/freq1[[3]]] // N

Fit Good distribution to

vest = 
      1/((freq1[[1]]/freq1[[2]])^(1/sest) - 1)

zest = ((vest + 3)/(vest + 2))^sest freq1[[4]]/
          freq1[[3]]

Method of moments

FitUncorrectedMoments[LerchDistribution[zest, sest, vest], freq1, 
      Frequency -> True, Categories -> classes1, 
      MaxIterations -> 100]

fitmom = 
      FitUncorrectedMoments[LerchDistribution[zest, sest, vest], data, 
        Frequency -> False, MaxIterations -> 100]

MMVarianceCovariance[
          LerchDistribution[0.7410652358091605`, 0.6254145541503756`, 
            0.6250698027028831`]]/1000 // MatrixForm

momfitfreq = 
      Map[PDF[Apply[LerchDistribution, fitmom], #] &, classes1] Length[
          data]

momfitplot = 
        ListPlot[Transpose[{classes1, momfitfreq}], PlotJoined -> True, 
          PlotStyle -> {Hue[0.65]}, 
          DisplayFunction -> Identity];

Show[Histogram[data, DisplayFunction -> Identity], momfitplot, 
        DisplayFunction -> $DisplayFunction];

putfreq[freq1, 
      classes1, {0, 1, 2, 3, 4, 5, 6, 7, 8, 10, 12, 14, 17}]

PearsonChiSquaredStatistic[Apply[LerchDistribution, fitmom], 
      freq1, {0, 1, 2, 3, 4, 5, 6, 7, 8, 10, 12, 14, 17}, 
      Frequency -> True, Categories -> classes1]

1 - 
      CDF[ChiSquareDistribution[12 - 3 - 1], 17.37190122707456`]

Maximum likelihood method

FitMLE[LerchDistribution[zest, sest, vest], data, 
      Frequency -> False, MaxIterations -> 100]

fitmle = 
      FitMLE[LerchDistribution[zest, sest, vest], freq1, 
        Frequency -> True, Categories -> classes1, 
        MaxIterations -> 100]

MLEVarianceCovariance[
          LerchDistribution[0.7356594726707477`, 0.5901541814054992`, 
            0.587518070064451`]]/1000 // MatrixForm

mlefitfreq = 
      Map[PDF[Apply[LerchDistribution, fitmle], #] &, classes1] Length[
          data]

mlefitplot = 
        ListPlot[Transpose[{classes1, mlefitfreq}], PlotJoined -> True, 
          PlotStyle -> {Hue[0.65]}, 
          DisplayFunction -> Identity];

Show[datahist, mlefitplot, 
        DisplayFunction -> $DisplayFunction];

putfreq[freq1, 
      classes1, {0, 1, 2, 3, 4, 5, 6, 7, 8, 10, 12, 14, 17}]

PearsonChiSquaredStatistic[Apply[LerchDistribution, fitmle], 
      freq1, {0, 1, 2, 3, 4, 5, 6, 7, 8, 10, 12, 14, 17}, 
      Frequency -> True, Categories -> classes1]

1 - 
      CDF[ChiSquareDistribution[12 - 3 - 1], 17.709318512547892`]

Example 5 (Minimize Chi-square)

Exmaple 5 shows how to fit data with Lerch distribution using minimization of chi^2.
Note: To fit special cases, Zipf, Zipf-Mandelbrot, and Good, specify 1/0 where appropriate as initial settings for parameters..

NOTE: please generate data in Example 3  first!!

Fit data with pearson's chi squared statistic...

putfreq[freq1, 
      classes1, {0, 1, 2, 3, 4, 5, 6, 7, 8, 11, 14, 17}]

FitPearsonChiSquared[LerchDistribution[0.7, 0.5, 0.5], 
      freq1, {0, 1, 2, 3, 4, 5, 6, 7, 8, 11, 14, 17}, Frequency -> True, 
      Categories -> classes1]

FitPearsonChiSquared[LerchDistribution[0.7, 0.5, 0.5], 
      data, {0, 1, 2, 3, 4, 5, 6, 7, 8, 11, 14, 17}]

PearsonChiSquaredStatistic[
      LerchDistribution[0.7841240834918685`, 0.8923115854901538`, 
        0.9552299791650437`], freq1, {0, 1, 2, 3, 4, 5, 6, 7, 8, 11, 14, 17}, 
      Frequency -> True, Categories -> classes1]

1 - 
      CDF[ChiSquareDistribution[12 - 3 - 1], 11.658881046492375`]

Options

Frequency

?Frequency

Categories

?Categories

MethodPhi

?MethodPhi