(*****************************************************************************
Conditions of Use

Of the Mathematica package LerchDistribution.m

by Sergej V. Aksenov (http://aksenov.freeshell.org/)


This software (file LerchDistribution.m) is copyright by Sergej V. Aksenov 
and Michael A. Savageau.

This software and its constituent parts come with no warranty, whether expressed 
or implied, that it is free of errors or suitable for any specific purpose.  It 
must not be used to solve any problem, whose incorrect solution could result in 
injury to a person, institution, or property.  It is user's own risk to use this 
software or its parts and the authors disclaim all liability for such use.

This software is distributed ``as is''. In particular, no maintenance, support, 
troubleshooting, or subsequent upgrade is implied.

The use of this software must be acknowledged in any publication that contains 
results obtained with it; please contact the author for the citation.

The free use of this software and its parts is restricted for research 
purposes; commercial use require permission and licensing from Sergej V. 
Aksenov and Michael A. Savageau.
******************************************************************************)



(*:Mathematica Version: 4.1 *)

(*:Package Version: 2.0 *)

(*:Name: LerchDistribution` *)

(*:Context: LerchDistribution`LerchDistribution` *)

(*:Title: Discrete Lerch distribution *)

(*:Author: Sergej V. Aksenov (http://aksenov.freeshell.org/) *)

(*:History:
   Created: November, 2001.
   v1.0; Last modified: November, 2001
   v2.0; Last modified: June, 2002 (functions are extended to handle special cases Zipf, Zipf-Mandelbrot, Good, and truncated distributions; there is option to use CNCT algorithm for LerchPhi; added function to calculate value Pearson's chi^2 statitic and function to minimize it). *)

(*:Copyright: Copyright 2001-2002, Sergej V. Aksenov and Michael A. Savageau *)

(*:Reference: Usage messages *)

(*:Summary:
This package provides properties and functionals of the discrete Lerch distribution. *)

(*:Keywords: Discrete distribution, Lerch distribution, Lerch's transcendent *)

(*:Limitations: *)

(*:Note: 
To Do:
*  all functions for lerch, Zipf, Zipf-Mandelbrot, and Good, and their truncated versions
*  FitPearson.. and PearsonChi.. should accept any PDF.
*)

(*:Sources:
"Parameter estimation and random number generation from a Zipf-related Lerch distribution", by S.V. Aksenov and M.A. Savageau, Computational Statistics and Data Analysis, submitted, 2001; preprint MPS: Statistics/0206008 (http://www.mathpreprints.com/)
*)



BeginPackage["LerchDistribution`LerchDistribution`","Statistics`DescriptiveStatistics`","Statistics`DataManipulation`","Statistics`Common`DistributionsCommon`","Utilities`FilterOptions`"]
        
             
             
(* Extend the usage of Random.  Look for the indicated phrase
        to determine if this has already been done. *)

If[StringPosition[Random::usage, "specified statistical distribution"] === {},
Random::usage = Random::usage <> " " <>
"Random[distribution] gives a random number with the specified
statistical distribution."
]



(* Exported symbols *)

LerchDistribution::usage = 
"LerchDistribution[z, s, v, interval] represents Lerch distribution with parameters z, s, and v that is defined on an interval Interval[{a,b}].  If interval is omitted, default Interval[{0,Infinity}] is assumed.  Setting z=1 and v=0 specifies Zipf distribution (domain at least Interval[{1,Infinity}]);  z=1 specifies Zipf-Mandelbrot distribution;  v=0 specifies Good distribution (domain at least Interval[{1,Infinity}])."

DescendingFactorialMoment::usage = "DescendingFactorialMoment[distribution, r] gives r-th descending factorial moment of the specified statistical distribution."

UncorrectedMoment::usage = "UncorrectedMoment[distribution, r] gives r-th uncorrected moment of the specified statistical distribution."

FitUncorrectedMoments::usage = "FitUncorrectedMoments[distribution, data] gives a list of estimates of parameters of the statistical distribution, whose specified parameters represent an initial guess, that fits data, as determined from solution of equations equating sample with population uncorrected moments.  Data can be given as a list of raw data (turn Frequency->False, default) or as a list of empirical frequencies of distinct points (turn Frequency->True), in which case one should also supply to option Categories a list of [distinct] points corresponding to the list of frequencies."

FitMLE::usage = "FitMLE[distribution, data] gives a list of estimates of parameters of the statistical distribution, whose specified parameters represent an initial guess, that fits data by maximizing the likelihood function.  Data can be given as a list of raw data (turn Frequency->False, default) or as a list of empirical frequencies of distinct points (turn Frequency->True), in which case one should also supply to option Categories a list of [distinct] points corresponding to the list of frequencies."

FitPearsonChiSquared::usage="FitPearsonChiSquared[distribution, data, cutoffs] gives a list of estimates of parameters of the statistical distribution, whose specified parameters represent an initial guess, that fits data by minimizing Pearson's Chi squared statistic over cells specified by a list of cutoffs.    Data can be given as a list of raw data (turn Frequency->False, default) or as a list of empirical frequencies of distinct points (turn Frequency->True), in which case one should also supply to option Categories a list of [distinct] points corresponding to the list of frequencies."

PearsonChiSquaredStatistic::usage="FitPearsonChiSquared[distribution, data, cutoffs] gives value of Pearson's Chi squared statistic calculated with the specified statistical distribution over cells of data specified by a list of cutoffs.    Data can be given as a list of raw data (turn Frequency->False, default) or as a list of empirical frequencies of distinct points (turn Frequency->True), in which case one should also supply to option Categories a list of [distinct] points corresponding to the list of frequencies."

MLEVarianceCovariance::usage = "MLEVarianceCovariance[distribution] gives a variance-covariance matrix for MLE estimates of parameters of the specified statistical distribution that has parameters fitted to data with MLE."

MMVarianceCovariance::usage = "MMVarianceCovariance[distribution] gives a variance-covariance matrix for method of moments estimates of parameters of the specified statistical distribution that has parameters fitted to data with MM."


Frequency::usage = "Frequency is option for FitUncorrectedMoments, FitMLE, FitPearsonChiSquared, PearsonChiSquaredStatistic that specifies whether data is given as a list of empirical frequences of distinct points (True) in which case one should also supply to option Categories a list of [distinct] points corresponding to the list of frequencies, or raw data (False)."

Categories::usage = "Categories is option for FitUncorrectedMoments, FitMLE, FitPearsonChiSquared, PearsonChiSquaredStatistic that specifies the list of [distinct] points corresponding to the list of frequencies suppled as data, if Frequency->True"
 
MethodPhi::usage =
"Option MethodPhi specifies either LerchPhi or LerchPhiCNCT as algorithms for calculating Lerch's transcendent Phi."

AdjustedGeneralizedPoissonDistribution::usage="AdjustedGeneralizedPoissonDistribution[theta,alpha,phi] represents adjusted generalized Poisson distribution."

GeneralizedPoissonDistribution::usage="GeneralizedPoissonDistribution[theta,lambda] represents generalized Poisson distribution."


(* Extend usage messages for functions defined in DescriptiveStatistics.m.
  Look for the indicated phrase to determine if this has already been done. *)

If[StringPosition[Mean::usage, "specified statistical distribution"] === {},

ExpectedValue::usage = ExpectedValue::usage <> " " <>
"ExpectedValue[f, distribution] gives the expected value of the pure function
f with respect to the specified statistical distribution.
ExpectedValue[f, distribution, x] gives the expected value of the
function f of x with respect to the specified univariate statistical
distribution for x.";

Mean::usage = Mean::usage <> " " <>
"Mean[distribution] gives the mean of the specified statistical distribution.";

CentralMoment::usage = CentralMoment::usage <> " " <>
"CentralMoment[distribution, r] gives r-th central moment of the specified statistical distribution with respect to its mean.";
 
Variance::usage = Variance::usage <> " " <>
"Variance[distribution] gives the variance of the specified statistical
distribution.";
 
StandardDeviation::usage = StandardDeviation::usage <> " " <>
"StandardDeviation[distribution] gives the standard deviation of the
specified statistical distribution.";
 
Skewness::usage = Skewness::usage <> " " <>
"Skewness[distribution] gives the coefficient of skewness of the
specified statistical distribution.";
 
Kurtosis::usage = Kurtosis::usage <> " " <>
"Kurtosis[distribution] gives the coefficient of kurtosis of the
specified statistical distribution.";
 
Quantile::usage = Quantile::usage <> " " <>
"Quantile[distribution, q] gives the q-th quantile of the specified
statistical distribution."
]


prt=Unprotect[ParameterQ,Domain,PDF,CDF,Mean,Variance,StandardDeviation,CentralMoment,Skewness,Kurtosis,CharacteristicFunction,Quantile,Random,RandomArray]


If[Options[ParameterQ]=!={}, Options[ParameterQ] = Join[Options[ParameterQ],{MethodPhi -> LerchPhi}], Options[ParameterQ] = {MethodPhi -> LerchPhi}]
If[Options[Domain]=!={}, Options[Domain] = Join[Options[Domain],{MethodPhi -> LerchPhi}], Options[Domain] = {MethodPhi -> LerchPhi}]
If[Options[PDF]=!={}, Options[PDF] = Join[Options[PDF],{MethodPhi -> LerchPhi}], Options[PDF] = {MethodPhi -> LerchPhi}]
If[Options[CDF]=!={}, Options[CDF] = Join[Options[CDF],{MethodPhi -> LerchPhi}], Options[CDF] = {MethodPhi -> LerchPhi}]
If[Options[Mean]=!={}, Options[Mean] = Join[Options[Mean],{MethodPhi -> LerchPhi}], Options[Mean] = {MethodPhi -> LerchPhi}]
If[Options[Variance]=!={}, Options[Variance] = Join[Options[Variance],{MethodPhi -> LerchPhi}], Options[Variance] = {MethodPhi -> LerchPhi}]
If[Options[StandardDeviation]=!={}, Options[StandardDeviation] = Join[Options[StandardDeviation],{MethodPhi -> LerchPhi}], Options[StandardDeviation] = {MethodPhi -> LerchPhi}]
If[Options[UncorrectedMoment]=!={}, Options[UncorrectedMoment] = Join[Options[UncorrectedMoment],{MethodPhi -> LerchPhi}], Options[UncorrectedMoment] = {MethodPhi -> LerchPhi}]
If[Options[ExpectedValue]=!={}, Options[ExpectedValue] = Join[Options[ExpectedValue],{MethodPhi -> LerchPhi}], Options[ExpectedValue] = {MethodPhi -> LerchPhi}]
If[Options[CentralMoment]=!={}, Options[CentralMoment] = Join[Options[CentralMoment],{MethodPhi -> LerchPhi}], Options[CentralMoment] = {MethodPhi -> LerchPhi}]
If[Options[DescendingFactorialMoment]=!={}, Options[DescendingFactorialMoment] = Join[Options[DescendingFactorialMoment],{MethodPhi -> LerchPhi}], Options[DescendingFactorialMoment] = {MethodPhi -> LerchPhi}]
If[Options[Skewness]=!={}, Options[Skewness] = Join[Options[Skewness],{MethodPhi -> LerchPhi}], Options[Skewness] = {MethodPhi -> LerchPhi}]
If[Options[Kurtosis]=!={}, Options[Kurtosis] = Join[Options[Kurtosis],{MethodPhi -> LerchPhi}], Options[Kurtosis] = {MethodPhi -> LerchPhi}]
If[Options[CharacteristicFunction]=!={}, Options[CharacteristicFunction] = Join[Options[CharacteristicFunction],{MethodPhi -> LerchPhi}], Options[CharacteristicFunction] = {MethodPhi -> LerchPhi}]
If[Options[Quantile]=!={}, Options[Quantile] = Join[Options[Quantile],{MethodPhi -> LerchPhi}], Options[Quantile] = {MethodPhi -> LerchPhi}]
If[Options[Random]=!={}, Options[Random] = Join[Options[Random],{MethodPhi -> LerchPhi}], Options[Random] = {MethodPhi -> LerchPhi}]
If[Options[RandomArray]=!={}, Options[RandomArray] = Join[Options[RandomArray],{MethodPhi -> LerchPhi}], Options[RandomArray] = {MethodPhi -> LerchPhi}]
If[Options[MLEVarianceCovariance]=!={}, Options[MLEVarianceCovariance] = Join[Options[MLEVarianceCovariance],{MethodPhi -> LerchPhi}], Options[MLEVarianceCovariance] = {MethodPhi -> LerchPhi}]
If[Options[MMVarianceCovariance]=!={}, Options[MMVarianceCovariance] = Join[Options[MMVarianceCovariance],{MethodPhi -> LerchPhi}], Options[MMVarianceCovariance] = {MethodPhi -> LerchPhi}]
If[Options[FitUncorrectedMoments]=!={}, Options[FitUncorrectedMoments] = Join[Options[FitUncorrectedMoments],{Frequency -> False, Categories -> {}, MethodPhi -> LerchPhi}], Options[FitUncorrectedMoments] = {Frequency -> False, Categories -> {}, MethodPhi -> LerchPhi}]
If[Options[FitMLE]=!={}, Options[FitMLE] = Join[Options[FitMLE],{Frequency -> False, Categories -> {}, MethodPhi -> LerchPhi}], Options[FitMLE] = {Frequency -> False, Categories -> {}, MethodPhi -> LerchPhi}]
If[Options[FitPearsonChiSquared]=!={}, Options[FitPearsonChiSquared] = Join[Options[FitPearsonChiSquared],{Frequency -> False, Categories -> {}, MethodPhi -> LerchPhi}], Options[FitPearsonChiSquared] = {Frequency -> False, Categories -> {}, MethodPhi -> LerchPhi}]
If[Options[PearsonChiSquaredStatistic]=!={}, Options[PearsonChiSquaredStatistic] = Join[Options[PearsonChiSquaredStatistic],{Frequency -> False, Categories -> {}, MethodPhi -> LerchPhi}], Options[PearsonChiSquaredStatistic] = {Frequency -> False, Categories -> {}, MethodPhi -> LerchPhi}] 



Begin["`Private`"]



Default[LerchDistribution,4] = Interval[{0,Infinity}]



LerchDistribution/: ParameterQ[LerchDistribution[z_,s_,v_,int_.], opts___?OptionQ] := And[
        If[FreeQ[N[z], Complex], True,
           Message[LerchDistribution::realparm, z]; False],
        If[FreeQ[N[s], Complex], True,
           Message[LerchDistribution::realparm, s]; False],
        If[FreeQ[N[v], Complex], True,
           Message[LerchDistribution::realparm, v]; False],
        If[
        (* Lerch *)
        (0<z<1 && (Min[int]+v)>0) || 
        (* Zipf-Mandelbrot [generalized zeta] *)
        (z==1 && v>0 && s>1) ||
        (* Zipf [zeta] *) 
        (z==1 && v==0 && Min[int]>=1 && s>1) ||
        (* Good [Jonquiere, aka PolyLog] *)
        (0<z<1 && v==0 && Min[int]>=1),
        True, Message[LerchDistribution::defpar,z,s,v]; False], True] /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && (MethodPhi/.{opts}/.Options[ParameterQ])===LerchPhi


LerchDistribution/: ParameterQ[LerchDistribution[z_,s_,v_,int_.], opts___?OptionQ] := And[
        If[FreeQ[N[z], Complex], True,
           Message[LerchDistribution::realparm, z]; False],
        If[FreeQ[N[s], Complex], True,
           Message[LerchDistribution::realparm, s]; False],
        If[FreeQ[N[v], Complex], True,
           Message[LerchDistribution::realparm, v]; False],
        If[
        (* Lerch only*)
        (0<z<1 && v>0),
        True, Message[LerchDistribution::defpar,z,s,v]; False], True] /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && (MethodPhi/.{opts}/.Options[ParameterQ])===Global`LerchPhiCNCT

LerchDistribution::realparm = "Parameter `` is expected to be real."
LerchDistribution::defpar = "LerchDistribution is not defined for z=`1`, s=`2`, v=`3`"


        
LerchDistribution/: DomainQ[LerchDistribution[z_,s_,v_,int_.], list_?VectorQ] :=
	Scan[If[!IntervalMemberQ[int,#], Return[False]]&, list] =!= False /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]]
	
LerchDistribution/: DomainQ[LerchDistribution[z_,s_,v_,int_.], x_] :=
	IntervalMemberQ[int,x] /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]]

LerchDistribution/: Domain[LerchDistribution[z_,s_,v_,int_.],opts___?OptionQ]:=int /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z,s,v,int],opts]



(* PDF *)

(* out of range or nonintegral *)

LerchDistribution/: PDF[LerchDistribution[z_,s_,v_,int_.], x_, opts___?OptionQ] := 0 /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && NumberQ[N[x]] && !DomainQ[LerchDistribution[z,s,v,int],x]

(* covers Lerch and Good *)

LerchDistribution/: PDF[LerchDistribution[z_,s_,v_,int_.], x_, opts___?OptionQ] := Module[{c,a,b}, a=Min[int]; b=Max[int]; If[b==Infinity, c=z^a LerchPhi[z,s,v+a], c=z^a LerchPhi[z,s,v+a] - z^(b+1) LerchPhi[z,s,v+b+1]];
   1.0/c z^x/(v+x)^s ]/; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z,s,v,int],opts] && Head[x] =!= List && (MethodPhi/.{opts}/.Options[PDF])===LerchPhi && z!=1 

(* Weird "bug" in the above function: inserting checking for x to be in the domain [using domainq] results in incorrect calculation in Sum (not 1 for z close to 1) *)

(* covers Zipf and Zipf-Mandelbrot *)

LerchDistribution/: PDF[LerchDistribution[z_,s_,v_,int_.], x_, opts___?OptionQ] := Module[{c,a,b}, a=Min[int]; b=Max[int]; If[b==Infinity, c=LerchPhi[1,s,v+a], c=LerchPhi[1,s,v+a]-LerchPhi[1,s,v+b+1]]; 
   1/(v+x)^s/c ]/; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z,s,v,int],opts] && Head[x] =!= List && (MethodPhi/.{opts}/.Options[PDF])===LerchPhi && z==1

(* covers Lerch only [with CNCT for Phi] *)

LerchDistribution/: PDF[LerchDistribution[z_,s_,v_,int_.], x_, opts___?OptionQ] := Module[{c,a,b}, a=Min[int]; b=Max[int]; If[b==Infinity, c=z^a Global`LerchPhiCNCT[z,s,v+a,opts], c=z^a Global`LerchPhiCNCT[z,s,v+a,opts] - z^(b+1) Global`LerchPhiCNCT[z,s,v+b+1,opts]]; 
   z^x/(v+x)^s/c ]/; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && (MethodPhi/.{opts}/.Options[PDF])===Global`LerchPhiCNCT && ParameterQ[LerchDistribution[z,s,v,int],opts] && Head[x] =!= List
        
        
        
(* CDF *)

(* zero sum *)

LerchDistribution/: CDF[LerchDistribution[z_,s_,v_,int_.], x_, opts___?OptionQ] := Which[x<Min[int],0,x>Max[int],1] /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z,s,v,int],opts] && NumberQ[x] && !DomainQ[LerchDistribution[z,s,v,int],x]

(* Covers Lerch and Good *)

LerchDistribution/: CDF[LerchDistribution[z_,s_,v_,int_.], x_, opts___?OptionQ] := Module[{a,b,c}, a = Min[int]; b = Max[int]; If[b==Infinity, c=z^a LerchPhi[z,s,v+a], c=z^a LerchPhi[z,s,v+a] - z^(b+1) LerchPhi[z,s,v+b+1]]; If[x==Infinity, 1, (z^a LerchPhi[z,s,v+a] - z^(x+1) LerchPhi[z,s,v+x+1]) / c]] /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z,s,v,int],opts] && Head[x] =!= List && (MethodPhi/.{opts}/.Options[CDF])===LerchPhi && z!=1

(* Covers Zipf and Zipf-Mandelbrot *)

LerchDistribution/: CDF[LerchDistribution[z_,s_,v_,int_.], x_, opts___?OptionQ] := Module[{a,b,c}, a = Min[int]; b = Max[int]; If[b==Infinity, c=LerchPhi[1,s,v+a], c=LerchPhi[1,s,v+a] - LerchPhi[1,s,v+b+1]]; If[x==Infinity, 1, (LerchPhi[1,s,v+a] - LerchPhi[1,s,v+x+1]) / c]] /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z,s,v,int],opts] && Head[x] =!= List && (MethodPhi/.{opts}/.Options[CDF])===LerchPhi && z==1

(* Covers Lerch only (via CNCT) *)

LerchDistribution/: CDF[LerchDistribution[z_,s_,v_,int_.], x_, opts___?OptionQ] := Module[{a,b,c}, a = Min[int]; b = Max[int]; If[b==Infinity, c=z^a Global`LerchPhiCNCT[z,s,v+a,opts], c=z^a Global`LerchPhiCNCT[z,s,v+a,opts] - z^(b+1) Global`LerchPhiCNCT[z,s,v+b+1,opts]];
If[x==Infinity, 1, (z^a Global`LerchPhiCNCT[z,s,v+a,opts] - z^(x+1) Global`LerchPhiCNCT[z,s,v+x+1,opts]) / c]] /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z,s,v,int],opts] && Head[x] =!= List && IntervalMemberQ[int,x] && (MethodPhi/.{opts}/.Options[CDF])===Global`LerchPhiCNCT



LerchDistribution::exmom="Moment of the Zipf or Zipf-Mandelbrot distribution for s=`` does not exist."


        
(* Mean *)

(* Covers Lerch and Good *)

LerchDistribution/: Mean[LerchDistribution[z_,s_,v_,int_.], opts___?OptionQ] := 
Module[{a,b,c,bterm1,bterm2}, a = Min[int]; b = Max[int]; If[b==Infinity, c=z^a LerchPhi[z,s,v+a]; bterm1=bterm2=0, c=z^a LerchPhi[z,s,v+a] - z^(b+1) LerchPhi[z,s,v+b+1]; bterm1=z^(b+1) LerchPhi[z,s-1,v+b+1]; bterm2=v z^(b+1) LerchPhi[z,s,v+b+1]];
(z^a LerchPhi[z,s-1,v+a]-bterm1-v z^a LerchPhi[z,s,v+a]+bterm2)/c
]/; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z,s,v,int],opts] && (MethodPhi/.{opts}/.Options[Mean])===LerchPhi && z!=1

(* Covers Zipf and Zipf-Mandelbrot *)

LerchDistribution/: Mean[LerchDistribution[z_,s_,v_,int_.], opts___?OptionQ] := 
Module[{a,b,c,bterm1,bterm2}, If[s<=2, Message[LerchDistribution::exmom,s]; Return[$Failed]]; a = Min[int]; b = Max[int]; If[b==Infinity, c=LerchPhi[1,s,v+a]; bterm1=bterm2=0, c=LerchPhi[1,s,v+a] - LerchPhi[1,s,v+b+1]; bterm1=LerchPhi[1,s-1,v+b+1]; bterm2=v LerchPhi[1,s,v+b+1]];
(LerchPhi[1,s-1,v+a]-bterm1-v LerchPhi[1,s,v+a]+bterm2)/c
]/; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z,s,v,int],opts] && (MethodPhi/.{opts}/.Options[Mean])===LerchPhi && z==1

(* Covers Lerch only (via CNCT) *)

LerchDistribution/: Mean[LerchDistribution[z_,s_,v_,int_.], opts___?OptionQ] := 
Module[{a,b,c,bterm1,bterm2}, a = Min[int]; b = Max[int]; If[b==Infinity, c=z^a Global`LerchPhiCNCT[z,s,v+a,opts]; bterm1=bterm2=0, c=z^a Global`LerchPhiCNCT[z,s,v+a,opts] - z^(b+1) Global`LerchPhiCNCT[z,s,v+b+1,opts]; bterm1=z^(b+1) Global`LerchPhiCNCT[z,s-1,v+b+1,opts]; bterm2=v z^(b+1) Global`LerchPhiCNCT[z,s,v+b+1,opts]];
(z^a Global`LerchPhiCNCT[z,s-1,v+a,opts]-bterm1-v z^a Global`LerchPhiCNCT[z,s,v+a,opts]+bterm2)/c
]/; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z,s,v,int],opts] && (MethodPhi/.{opts}/.Options[Mean])===Global`LerchPhiCNCT



(* UncorrectedMoment *)

(* Covers Lerch and Good *)

LerchDistribution/: UncorrectedMoment[LerchDistribution[z_,s_,v_,int_.],r_, opts___?OptionQ] := Module[{j,a,b,c},
a = Min[int]; b = Max[int]; If[b==Infinity, c=z^a LerchPhi[z,s,v+a], c=z^a LerchPhi[z,s,v+a] - z^(b+1) LerchPhi[z,s,v+b+1]];
If[v==0,(z^a LerchPhi[z,s-r,a]-z^(b+1) LerchPhi[z,s-r,b+1])/c,Sum[ Binomial[r, j] (-v)^(r-j) (z^a LerchPhi[z,s-j,v+a]-z^(b+1) LerchPhi[z,s-j,v+b+1]), {j,0,r}] / c]] /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z,s,v,int],opts] && (MethodPhi/.{opts}/.Options[UncorrectedMoment])===LerchPhi && z!=1

(* Covers Zipf and Zipf-Mandelbrot *)

LerchDistribution/: UncorrectedMoment[LerchDistribution[z_,s_,v_,int_.],r_, opts___?OptionQ] := Module[{j,a,b,c,bterm1,bterm2}, If[s<=r+1, Message[LerchDistribution::exmom,s]; Return[$Failed]];
a = Min[int]; b = Max[int]; If[b==Infinity, c=LerchPhi[1,s,v+a]; bterm1=bterm2=0, c=LerchPhi[1,s,v+a] - LerchPhi[1,s,v+b+1]; bterm1=LerchPhi[1,s-r,b+1]; bterm2=LerchPhi[1,s-j,v+b+1]];
If[v==0,(LerchPhi[1,s-r,a]-bterm1)/c,Sum[ Binomial[r, j] (-v)^(r-j) (LerchPhi[1,s-j,v+a]-bterm2), {j,0,r}] / c]] /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z,s,v,int],opts] && (MethodPhi/.{opts}/.Options[UncorrectedMoment])===LerchPhi && z==1

(* Covers Lerch only (via CNCT) *)

LerchDistribution/: UncorrectedMoment[LerchDistribution[z_,s_,v_,int_.],r_, opts___?OptionQ] := Module[{j,a,b,c},
a = Min[int]; b = Max[int]; If[b==Infinity, c=z^a Global`LerchPhiCNCT[z,s,v+a,opts], c=z^a Global`LerchPhiCNCT[z,s,v+a,opts] - z^(b+1) Global`LerchPhiCNCT[z,s,v+b+1,opts]];
If[v==0,(z^a Global`LerchPhiCNCT[z,s-r,a,opts]-z^(b+1) Global`LerchPhiCNCT[z,s-r,b+1,opts])/c,Sum[ Binomial[r, j] (-v)^(r-j) (z^a Global`LerchPhiCNCT[z,s-j,v+a,opts]-z^(b+1) Global`LerchPhiCNCT[z,s-j,v+b+1,opts]), {j,0,r}] / c]] /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z,s,v,int],opts] && (MethodPhi/.{opts}/.Options[UncorrectedMoment])===Global`LerchPhiCNCT





LerchDistribution/: Skewness[LerchDistribution[z_,s_,v_,int_.], opts___?OptionQ] := Module[{t0,t1,t2,t3}, t0 = LerchPhi[z,s,v]; t1 = LerchPhi[z,s-1,v]; t2 = LerchPhi[z,s-2,v]; t3 = LerchPhi[z,s-3,v]; mu3 = -(t1/t0)^3 + ( 3 t1^3/t0^2 - 3 t1 t2/t0 + t3)/t0; mu2 = (t1/t0)^2 + ( -2 t1^2/t0 + t2)/t0; Return[mu3/mu2^1.5]] /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z,s,v,int],opts] && (MethodPhi/.{opts}/.Options[Skewness])===LerchPhi

LerchDistribution/: Kurtosis[LerchDistribution[z_,s_,v_,int_.], opts___?OptionQ] := Module[{t0,t1,t2,t3,t4}, t0 = LerchPhi[z,s,v]; t1 = LerchPhi[z,s-1,v]; t2 = LerchPhi[z,s-2,v]; t3 = LerchPhi[z,s-3,v]; t4 = LerchPhi[z,s-4,v]; mu2 = (t1/t0)^2 + ( -2 t1^2/t0 + t2)/t0; mu4 = (t1/t0)^4 + ( -4 t1^4/t0^3 + 6 (t1/t0)^2 t2 - 4 t1/t0 t3 + t4)/t0; Return[mu4/mu2^2]] /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z,s,v,int],opts] && (MethodPhi/.{opts}/.Options[Kurtosis])===LerchPhi






(* CentralMoment *)

(* Covers Lerch and Good *)

LerchDistribution/: CentralMoment[LerchDistribution[z_,s_,v_,int_.],r_, opts___?OptionQ] := Module[{j,a,b,c,c1,vpmu},
a = Min[int]; b = Max[int]; If[b==Infinity, c=z^a LerchPhi[z,s,v+a]; c1=z^a LerchPhi[z,s-1,v+a], c=z^a LerchPhi[z,s,v+a] - z^(b+1) LerchPhi[z,s,v+b+1]; c1=z^a LerchPhi[z,s-1,v+a] - z^(b+1) LerchPhi[z,s-1,v+b+1]]; vpmu=c1/c;
If[b==Infinity, Sum[ Binomial[r, j] (-vpmu)^(r-j) z^a LerchPhi[z,s-j,v+a], {j,0,r}] / c, Sum[ Binomial[r, j] (-vpmu)^(r-j) (z^a LerchPhi[z,s-j,v+a]-z^(b+1) LerchPhi[z,s-j,v+b+1]), {j,0,r}] / c]] /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z,s,v,int],opts] && (MethodPhi/.{opts}/.Options[CentralMoment])===LerchPhi && z!=1

(* Covers Zipf and Zipf-Mandelbrot *)

LerchDistribution/: CentralMoment[LerchDistribution[z_,s_,v_,int_.],r_, opts___?OptionQ] := Module[{j,a,b,c,c1,vpmu,bterm}, If[s<=r+1, Message[LerchDistribution::exmom,s]; Return[$Failed]];
a = Min[int]; b = Max[int]; If[b==Infinity, bterm=0; c=LerchPhi[1,s,v+a]; c1=LerchPhi[1,s-1,v+a], bterm=LerchPhi[z,s-j,v+b+1]; c=LerchPhi[1,s,v+a] - LerchPhi[1,s,v+b+1]; c1=LerchPhi[1,s-1,v+a] - LerchPhi[1,s-1,v+b+1]]; vpmu=c1/c;
Sum[ Binomial[r, j] (-vpmu)^(r-j) (LerchPhi[1,s-j,v+a]-bterm), {j,0,r}] / c] /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z,s,v,int],opts] && (MethodPhi/.{opts}/.Options[CentralMoment])===LerchPhi && z==1

(* Covers Lerch only (via CNCT) *)

LerchDistribution/: CentralMoment[LerchDistribution[z_,s_,v_,int_.],r_, opts___?OptionQ] := Module[{j,a,b,c,c1,vpmu},
a = Min[int]; b = Max[int]; If[b==Infinity, c=z^a Global`LerchPhiCNCT[z,s,v+a,opts]; c1=z^a Global`LerchPhiCNCT[z,s-1,v+a,opts], c=z^a Global`LerchPhiCNCT[z,s,v+a,opts] - z^(b+1) Global`LerchPhiCNCT[z,s,v+b+1,opts]; c1=z^a Global`LerchPhiCNCT[z,s-1,v+a,opts] - z^(b+1) Global`LerchPhiCNCT[z,s-1,v+b+1,opts]]; vpmu=c1/c;
Sum[ Binomial[r, j] (-vpmu)^(r-j) (z^a Global`LerchPhiCNCT[z,s-j,v+a,opts]-z^(b+1) Global`LerchPhiCNCT[z,s-j,v+b+1,opts]), {j,0,r}] / c] /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z,s,v,int],opts] && (MethodPhi/.{opts}/.Options[CentralMoment])===LerchPhiCNCT



(* Variance *)

(* Covers Lerch and Good *)

LerchDistribution/: Variance[LerchDistribution[z_,s_,v_,int_.], opts___?OptionQ] := Module[{j,a,b,c,c1,d1,d2},
a = Min[int]; b = Max[int]; If[b==Infinity, c=z^a LerchPhi[z,s,v+a]; c1=z^a LerchPhi[z,s-1,v+a]; d1=c1; d2=z^a LerchPhi[z,s-2,v+a], c=z^a LerchPhi[z,s,v+a] - z^(b+1) LerchPhi[z,s,v+b+1]; c1=z^a LerchPhi[z,s-1,v+a] - z^(b+1) LerchPhi[z,s-1,v+b+1]; d1=z^a LerchPhi[z,s-1,v+a]-z^(b+1) LerchPhi[z,s-1,v+b+1]; d2=z^a LerchPhi[z,s-2,v+a]-z^(b+1) LerchPhi[z,s-2,v+b+1]]; 
((c1^2-2 c1 d1)/c+d2) / c] /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z,s,v,int],opts] && (MethodPhi/.{opts}/.Options[Variance])===LerchPhi && z!=1

(* Covers Zipf and Zipf-Mandelbrot *)

LerchDistribution/: Variance[LerchDistribution[z_,s_,v_,int_.], opts___?OptionQ] := Module[{j,a,b,c,c1,d1,d2}, If[s<=3, Message[LerchDistribution::exmom,s]; Return[$Failed]];
a = Min[int]; b = Max[int]; If[b==Infinity, c=LerchPhi[1,s,v+a]; c1=LerchPhi[1,s-1,v+a]; d1=c1; d2=LerchPhi[1,s-2,v+a], c=LerchPhi[1,s,v+a] - LerchPhi[1,s,v+b+1]; c1=LerchPhi[1,s-1,v+a] - LerchPhi[1,s-1,v+b+1]; d1=LerchPhi[1,s-1,v+a]-LerchPhi[1,s-1,v+b+1]; d2=LerchPhi[1,s-2,v+a]-LerchPhi[1,s-2,v+b+1]]; 
((c1^2-2 c1 d1)/c+d2) / c] /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z,s,v,int],opts] && (MethodPhi/.{opts}/.Options[Variance])===LerchPhi && z==1

(* Covers Lerch only (via CNCT) *)

LerchDistribution/: Variance[LerchDistribution[z_,s_,v_,int_.], opts___?OptionQ] := Module[{j,a,b,c,c1,d1,d2},
a = Min[int]; b = Max[int]; If[b==Infinity, c=z^a Global`LerchPhiCNCT[z,s,v+a,opts]; c1=z^a Global`LerchPhiCNCT[z,s-1,v+a,opts]; d1=c1; d2=z^a Global`LerchPhiCNCT[z,s-2,v+a,opts], c=z^a Global`LerchPhiCNCT[z,s,v+a,opts] - z^(b+1) Global`LerchPhiCNCT[z,s,v+b+1,opts]; c1=z^a Global`LerchPhiCNCT[z,s-1,v+a,opts] - z^(b+1) Global`LerchPhiCNCT[z,s-1,v+b+1,opts]; d1=z^a Global`LerchPhiCNCT[z,s-1,v+a,opts]-z^(b+1) Global`LerchPhiCNCT[z,s-1,v+b+1,opts]; d2=z^a Global`LerchPhiCNCT[z,s-2,v+a,opts]-z^(b+1) Global`LerchPhiCNCT[z,s-2,v+b+1,opts]]; 
((c1^2-2 c1 d1)/c+d2) / c] /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z,s,v,int],opts] && (MethodPhi/.{opts}/.Options[Variance])===Global`LerchPhiCNCT
		
		
		
(* StandardDeviation *)

(* Covers Lerch and Good *)

LerchDistribution/: StandardDeviation[LerchDistribution[z_,s_,v_,int_.], opts___?OptionQ] := Module[{j,a,b,c,c1,d1,d2},
a = Min[int]; b = Max[int]; If[b==Infinity, c=z^a LerchPhi[z,s,v+a]; c1=z^a LerchPhi[z,s-1,v+a]; d1=c1; d2=z^a LerchPhi[z,s-2,v+a], c=z^a LerchPhi[z,s,v+a] - z^(b+1) LerchPhi[z,s,v+b+1]; c1=z^a LerchPhi[z,s-1,v+a] - z^(b+1) LerchPhi[z,s-1,v+b+1]; d1=z^a LerchPhi[z,s-1,v+a]-z^(b+1) LerchPhi[z,s-1,v+b+1]; d2=z^a LerchPhi[z,s-2,v+a]-z^(b+1) LerchPhi[z,s-2,v+b+1]]; 
Sqrt[((c1^2-2 c1 d1)/c+d2) / c]] /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z,s,v,int],opts] && (MethodPhi/.{opts}/.Options[StandardDeviation])===LerchPhi && z!=1

(* Covers Zipf and Zipf-Mandelbrot *)

LerchDistribution/: StandardDeviation[LerchDistribution[z_,s_,v_,int_.], opts___?OptionQ] := Module[{j,a,b,c,c1,d1,d2}, If[s<=3, Message[LerchDistribution::exmom,s]; Return[$Failed]];
a = Min[int]; b = Max[int]; If[b==Infinity, c=LerchPhi[1,s,v+a]; c1=LerchPhi[1,s-1,v+a]; d1=c1; d2=LerchPhi[1,s-2,v+a], c=LerchPhi[1,s,v+a] - LerchPhi[1,s,v+b+1]; c1=LerchPhi[1,s-1,v+a] - LerchPhi[1,s-1,v+b+1]; d1=LerchPhi[1,s-1,v+a]-LerchPhi[1,s-1,v+b+1]; d2=LerchPhi[1,s-2,v+a]-LerchPhi[1,s-2,v+b+1]]; 
Sqrt[((c1^2-2 c1 d1)/c+d2) / c]] /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z,s,v,int],opts] && (MethodPhi/.{opts}/.Options[StandardDeviation])===LerchPhi && z==1

(* Covers Lerch only (via CNCT) *)

LerchDistribution/: StandardDeviation[LerchDistribution[z_,s_,v_,int_.], opts___?OptionQ] := Module[{j,a,b,c,c1,d1,d2},
a = Min[int]; b = Max[int]; If[b==Infinity, c=z^a Global`LerchPhiCNCT[z,s,v+a,opts]; c1=z^a Global`LerchPhiCNCT[z,s-1,v+a,opts]; d1=c1; d2=z^a Global`LerchPhiCNCT[z,s-2,v+a,opts], c=z^a Global`LerchPhiCNCT[z,s,v+a,opts] - z^(b+1) Global`LerchPhiCNCT[z,s,v+b+1,opts]; c1=z^a Global`LerchPhiCNCT[z,s-1,v+a,opts] - z^(b+1) Global`LerchPhiCNCT[z,s-1,v+b+1,opts]; d1=z^a Global`LerchPhiCNCT[z,s-1,v+a,opts]-z^(b+1) Global`LerchPhiCNCT[z,s-1,v+b+1,opts]; d2=z^a Global`LerchPhiCNCT[z,s-2,v+a,opts]-z^(b+1) Global`LerchPhiCNCT[z,s-2,v+b+1,opts]]; 
Sqrt[((c1^2-2 c1 d1)/c+d2) / c]] /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z,s,v,int],opts] && (MethodPhi/.{opts}/.Options[StandardDeviation])===Global`LerchPhiCNCT



LerchDistribution/: ExpectedValue[f_Function, LerchDistribution[z_,s_,v_,int_.], opts___?OptionQ] := Module[{x}, NSum[ f[x] PDF[LerchDistribution[z,s,v,int], x, opts], {x, 0, Infinity}]]  /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z,s,v,int],opts] && (MethodPhi/.{opts}/.Options[ExpectedValue])===LerchPhi

LerchDistribution/: ExpectedValue[f_Function, LerchDistribution[z_,s_,v_,int_.], opts___?OptionQ] := Module[{x}, NSum[ f[x] PDF[LerchDistribution[z,s,v,int], x, opts], {x, 0, Infinity}]]  /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z,s,v,int],opts] && (MethodPhi/.{opts}/.Options[ExpectedValue])===Global`LerchPhiCNCT

LerchDistribution/: ExpectedValue[f_, LerchDistribution[z_,s_,v_,int_.], x_Symbol, opts___?OptionQ] := Module[{xx}, Sum[f PDF[LerchDistribution[z,s,v,int], x,opts], {x, 0, Infinity}]]/; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z,s,v,int],opts] && (MethodPhi/.{opts}/.Options[ExpectedValue])===LerchPhi

LerchDistribution/: ExpectedValue[f_, LerchDistribution[z_,s_,v_,int_.], x_Symbol, opts___?OptionQ] := Module[{xx}, Sum[f PDF[LerchDistribution[z,s,v,int], x,opts], {x, 0, Infinity}]]/; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z,s,v,int],opts] && (MethodPhi/.{opts}/.Options[ExpectedValue])===Global`LerchPhiCNCT



LerchDistribution/: DescendingFactorialMoment[LerchDistribution[z_,s_,v_,int_.],r_, opts___?OptionQ] := Module[{i}, Sum[ StirlingS1[r,i] UncorrectedMoment[LerchDistribution[z,s,v,int],i,opts], {i,0,r}]] /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z,s,v,int],opts] && (MethodPhi/.{opts}/.Options[DescendingFactorialMoment])===LerchPhi

LerchDistribution/: DescendingFactorialMoment[LerchDistribution[z_,s_,v_,int_.],r_, opts___?OptionQ] := Module[{i}, Sum[ StirlingS1[r,i] UncorrectedMoment[LerchDistribution[z,s,v,int],i,opts], {i,0,r}]] /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z,s,v,int],opts] && (MethodPhi/.{opts}/.Options[DescendingFactorialMoment])===Global`LerchPhiCNCT



LerchDistribution/: CharacteristicFunction[
        LerchDistribution[z_,s_,v_,int_.], t_, opts___?OptionQ] := LerchPhi[z E^(I t),s,v] / LerchPhi[z,s,v] /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z,s,v,int],opts] && (MethodPhi/.{opts}/.Options[CharacteristicFunction])===LerchPhi

LerchDistribution/: CharacteristicFunction[
        LerchDistribution[z_,s_,v_,int_.], t_, opts___?OptionQ] := Global`LerchPhiCNCT[z E^(I t),s,v,opts] / Global`LerchPhiCNCT[z,s,v,opts] /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z,s,v,int],opts] && (MethodPhi/.{opts}/.Options[CharacteristicFunction])===Global`LerchPhiCNCT
				
				
				
LerchDistribution/: Quantile[LerchDistribution[z_,s_,v_,int_.], q_, opts___?OptionQ] := Module[{l0,l1,x}, l0 = LerchPhi[z,s,v]; l1 = LerchPhi[z,s,v+1] / l0; x = Floor[Log[(1 - q)/l1]/Log[z]]-1; If[x<0,x=0]; While[q <= 1 - z^(x+1) LerchPhi[z,s,v+x+1]/l0, x--]; Return[x+1]] /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && QuantileQ[q] && ParameterQ[LerchDistribution[z,s,v,int],opts] && (MethodPhi/.{opts}/.Options[Quantile])===LerchPhi

LerchDistribution/: Quantile[LerchDistribution[z_,s_,v_,int_.], q_, opts___?OptionQ] := Module[{l0,l1,x}, l0 = Global`LerchPhiCNCT[z,s,v,opts]; l1 = Global`LerchPhiCNCT[z,s,v+1,opts] / l0; x = Floor[Log[(1 - q)/l1]/Log[z]]-1; If[x<0,x=0]; While[q <= 1 - z^(x+1) Global`LerchPhiCNCT[z,s,v+x+1,opts]/l0, x--]; Return[x+1]] /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && QuantileQ[q] && ParameterQ[LerchDistribution[z,s,v,int],opts] && (MethodPhi/.{opts}/.Options[Quantile])===Global`LerchPhiCNCT
				
(* Inversion by correction *)

lerchrnst[z_, s_, v_, int_] := Module[{u, x, l0, l1}, l0 = LerchPhi[z,s,v]; l1 = LerchPhi[z,s,v+1] / l0; u = Random[]; x = Floor[Log[(1 - u)/l1]/Log[z]]-1; If[x<0,x=0]; If[s>0,While[u <= 1 - z^(x+1) LerchPhi[z,s,v+x+1]/l0, x--]; Return[x+1], While[u > 1 - z^(x+1) LerchPhi[z,s,v+x+1]/l0, x++]; Return[x]]]

lerchrncnct[z_, s_, v_, int_, opts___?OptionQ] := Module[{u, x, l0, l1}, l0 = Global`LerchPhiCNCT[z,s,v,opts]; l1 = Global`LerchPhiCNCT[z,s,v+1,opts] / l0; u = Random[]; x = Floor[Log[(1 - u)/l1]/Log[z]]-1; If[x<0,x=0]; If[s>0,While[u <= 1 - z^(x+1) Global`LerchPhiCNCT[z,s,v+x+1,opts]/l0, x--]; Return[x+1], While[u > 1 - z^(x+1) Global`LerchPhiCNCT[z,s,v+x+1,opts]/l0, x++]; Return[x]]]

(* Inversion by sequential search *)

(* lerchrn[z_, s_, v_, int_] := Module[{u, x, l0, l1}, l0 = LerchPhi[z,s,v]; l1 = LerchPhi[z,s,v+1] / l0; u = Random[]; x = 0; While[u > 1 - z^(x+1) LerchPhi[z,s,v+x+1]/l0, x++]; Return[x]] *)

LerchDistribution/: Random[LerchDistribution[z_,s_,v_,int_.], opts___?OptionQ] := lerchrnst[z,s,v,int] /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && (MethodPhi/.{opts}/.Options[Random])===LerchPhi

LerchDistribution/: Random[LerchDistribution[z_,s_,v_,int_.], opts___?OptionQ] := lerchrncnct[z,s,v,int,opts] /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && (MethodPhi/.{opts}/.Options[Random])===Global`LerchPhiCNCT



LerchDistribution/: RandomArray[LerchDistribution[z_,s_,v_,int_.], dim_, opts___?OptionQ] := Module[{m, array}, m = If[VectorQ[dim], Apply[Times, dim], dim]; array = Table[lerchrnst[z,s,v,int], {m}]; If[VectorQ[dim] && Length[dim] > 1, Fold[Partition[#1, #2]&, array, Reverse[Drop[dim, 1]] ], array]] /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ((IntegerQ[dim] && dim > 0) || VectorQ[dim, (IntegerQ[#] && # > 0)&]) && (MethodPhi/.{opts}/.Options[RandomArray])===LerchPhi

LerchDistribution/: RandomArray[LerchDistribution[z_,s_,v_,int_.], dim_, opts___?OptionQ] := Module[{m, array}, m = If[VectorQ[dim], Apply[Times, dim], dim]; array = Table[lerchrncnct[z,s,v,int,opts], {m}]; If[VectorQ[dim] && Length[dim] > 1, Fold[Partition[#1, #2]&, array, Reverse[Drop[dim, 1]] ], array]] /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ((IntegerQ[dim] && dim > 0) || VectorQ[dim, (IntegerQ[#] && # > 0)&]) && (MethodPhi/.{opts}/.Options[RandomArray])===Global`LerchPhiCNCT








LerchDistribution/: FitUncorrectedMoments[LerchDistribution[z1_, s1_, v1_, int_.], data1_?VectorQ, opts___?OptionQ] := Module[{m1,m2,m3,freq,cat,n,z,s,v,ret,ini,a,b,data}, 

   n = Length[data1];
   ini = {z1, s1, v1};

   freq = Frequency /. {opts} /. Options[FitUncorrectedMoments];
   If[(freq =!= True) && (freq =!= False), Message[FitUncorrectedMoments::type]; freq = False];

   cat = Categories /. {opts} /. Options[FitUncorrectedMoments];
   If[freq && VectorQ[cat] && !(Length[cat]>0 && Apply[And,Map[IntegerQ,cat]] && Apply[And,Map[NonNegative,cat]]), Message[FitUncorrectedMoments::cat]; Return[$Failed]];

   If[freq && Length[cat]>0 && Length[cat]<n, Message[FitUncorrectedMoments::len]; Return[$Failed]];

   If[freq, data = data1/Apply[Plus,data1],data=data1];

   If[freq, m1 = Apply[Plus, data cat]; m2 = Apply[Plus, data (cat)^2]; m3 = Apply[Plus, data (cat)^3], m1 = Apply[Plus, data]/n; m2 = Apply[Plus, data^2]/n; m3 = Apply[Plus, data^3]/n];

   With[{optf = FilterOptions[FindRoot, Sequence@@Flatten[{opts}]]}, ret = {z,s,v}/.FindRoot[{m1 == LerchPhi[z,s-1,v]/LerchPhi[z,s,v]-v, m2 == v^2 -2 v LerchPhi[z,s-1,v]/LerchPhi[z,s,v] + LerchPhi[z,s-2,v]/LerchPhi[z,s,v], m3 == -v^3 + 3 v^2 LerchPhi[z,s-1,v]/LerchPhi[z,s,v] - 3 v LerchPhi[z,s-2,v]/LerchPhi[z,s,v] + LerchPhi[z,s-3,v]/LerchPhi[z,s,v]}, {z,{ini[[1]],ini[[1]] 0.5}}, {s,{ini[[2]],ini[[2]] 0.5}}, {v,{ini[[3]],ini[[3]] 0.5}}, optf]];

   Return[ret];] /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z1,s1,v1,int],opts] && (MethodPhi/.{opts}/.Options[FitUncorrectedMoments])===LerchPhi


LerchDistribution/: FitUncorrectedMoments[LerchDistribution[z1_, s1_, v1_, int_.], data1_?VectorQ, opts___?OptionQ] := Module[{m1,m2,m3,freq,cat,n,z,s,v,ret,ini,data}, 

   n = Length[data1];
   ini = {z1, s1, v1};

   freq = Frequency /. {opts} /. Options[FitUncorrectedMoments];
   If[(freq =!= True) && (freq =!= False), Message[FitUncorrectedMoments::type]; freq = False];

   cat = Categories /. {opts} /. Options[FitUncorrectedMoments];
   If[freq && VectorQ[cat] && !(Length[cat]>0 && Apply[And,Map[IntegerQ,cat]] && Apply[And,Map[NonNegative,cat]]), Message[FitUncorrectedMoments::cat]; Return[$Failed]];

   If[freq && Length[cat]>0 && Length[cat]<n, Message[FitUncorrectedMoments::len]; Return[$Failed]];

   If[freq, data = data1/Apply[Plus,data1],data=data1];

   If[freq, m1 = Apply[Plus, data cat]; m2 = Apply[Plus, data (cat)^2]; m3 = Apply[Plus, data (cat)^3], m1 = Apply[Plus, data]/n; m2 = Apply[Plus, data^2]/n; m3 = Apply[Plus, data^3]/n];

   With[{optf = FilterOptions[FindRoot, Sequence@@Flatten[{opts}]]}, ret = {z,s,v}/.FindRoot[{m1 == Global`LerchPhiCNCT[z,s-1,v,opts]/Global`LerchPhiCNCT[z,s,v,opts]-v, m2 == v^2 -2 v Global`LerchPhiCNCT[z,s-1,v,opts]/Global`LerchPhiCNCT[z,s,v,opts] + Global`LerchPhiCNCT[z,s-2,v,opts]/Global`LerchPhiCNCT[z,s,v,opts], m3 == -v^3 + 3 v^2 Global`LerchPhiCNCT[z,s-1,v,opts]/Global`LerchPhiCNCT[z,s,v,opts] - 3 v Global`LerchPhiCNCT[z,s-2,v,opts]/Global`LerchPhiCNCT[z,s,v,opts] + Global`LerchPhiCNCT[z,s-3,v,opts]/Global`LerchPhiCNCT[z,s,v,opts]}, {z,{ini[[1]],ini[[1]] 0.5}}, {s,{ini[[2]],ini[[2]] 0.5}}, {v,{ini[[3]],ini[[3]] 0.5}}, optf]];

   Return[ret];] /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z1,s1,v1,int],opts] && (MethodPhi/.{opts}/.Options[FitUncorrectedMoments])===Global`LerchPhiCNCT

FitUncorrectedMoments::type = "Warning: option Frequency is not specified, using False."
FitUncorrectedMoments::cat = "Option Categories is not a vector of nonnegative integers."
FitUncorrectedMoments::len = "Lengths of data and Categories are incompatible."





LerchDistribution/: FitMLE[LerchDistribution[z1_, s1_, v1_, int_.], data1_?VectorQ, opts___?OptionQ] := Module[{e1,e2,e3,freq,cat,n,ret,i,z,s,v,ini,a,b,data}, 

   n = Length[data1];
   ini = {z1, s1, v1};

   freq = Frequency /. {opts} /. Options[FitMLE];
   If[(freq =!= True) && (freq =!= False), Message[FitMLE::type]; freq = False];

   cat = Categories /. {opts} /. Options[FitMLE];
   If[freq && VectorQ[cat] && !(Length[cat]>0 && Apply[And,Map[IntegerQ,cat]] && Apply[And,Map[NonNegative,cat]]), Message[FitMLE::cat]; Return[$Failed]];

   If[freq && Length[cat]>0 && Length[cat]<n, Message[FitMLE::len]; Return[$Failed]];

   If[freq, data = data1/Apply[Plus,data1],data=data1];

   If[freq, e1 = Apply[Plus, data cat]; e2 = (Apply[Plus, data Log[#+cat]])&; e3 = (Apply[Plus, data/(#+cat)])&, e1 = Apply[Plus, data]/n; e2 = (Apply[Plus, Log[#+data]]/n)&; e3 = (Apply[Plus, 1/(#+data)]/n)&];

   a = Min[int];  b = Max[int];
   If[freq && (a>Min[cat] || b<Max[cat]), Message[FitMLE::catf]; Return[$Failed]];
   
   rhs1[z_,s_,v_] := z/(z^a LerchPhi[z,s,v+a]) (z^a (a/z LerchPhi[z,s,v+a]+LerchPhi[z,s-1,v+a+1]-v LerchPhi[z,s,v+a+1])) /; b==Infinity;

   rhs1[z_,s_,v_] := z/(z^a LerchPhi[z,s,v+a]-z^(b+1) LerchPhi[z,s,v+b+1]) (z^a (a/z LerchPhi[z,s,v+a]+LerchPhi[z,s-1,v+a+1]-v LerchPhi[z,s,v+a+1]) - z^(b+1) ((b+1)/z LerchPhi[z,s,v+b+1]+LerchPhi[z,s-1,v+b+2]-v LerchPhi[z,s,v+b+2])) /; b!=Infinity;

   rhs2[z_,s_,v_] := Sum[z^i Log[v+a+i]/(v+a+i)^s,{i,0,Infinity}]/ LerchPhi[z,s,v+a] /; b==Infinity;

   rhs2[z_,s_,v_] := (z^a Sum[z^i Log[v+a+i]/(v+a+i)^s,{i,0,Infinity}]-z^(b+1) Sum[z^i Log[v+b+i+1]/(v+b+i+1)^s,{i,0,Infinity}])/(z^a LerchPhi[z,s,v+a]-z^(b+1) LerchPhi[z,s,v+b+1]) /; b!=Infinity;

   rhs3[z_,s_,v_] := LerchPhi[z,s+1,v+a]/LerchPhi[z,s,v+a] /; b==Infinity;

   rhs3[z_,s_,v_] := (z^a LerchPhi[z,s+1,v+a]-z^(b+1)/s LerchPhi[z,s+1,v+b+1])/(z^a LerchPhi[z,s,v+a]-z^(b+1) LerchPhi[z,s,v+b+1]) /; b!=Infinity;

(*   rhs1 = z (LerchPhi[z,s-1,v+1] - v LerchPhi[z,s,v+1])/LerchPhi[z,s,v];
   rhs2 = Sum[z^i Log[v+i] / (v+i)^s, {i,0,Infinity}]/LerchPhi[z,s,v];
   rhs3 = LerchPhi[z,s+1,v]/LerchPhi[z,s,v]; *)

   With[{optf = FilterOptions[FindRoot, Sequence@@Flatten[{opts}]]}, ret = {z,s,v}/.FindRoot[{e1 == rhs1[z,s,v], e2[v] == rhs2[z,s,v], e3[v] == rhs3[z,s,v]}, {z,{ini[[1]],ini[[1]] 0.9}}, {s,{ini[[2]],ini[[2]] 0.9}}, {v,{ini[[3]],ini[[3]] 0.9}}, optf]];

   Return[ret];] /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z1,s1,v1,int],opts] && (MethodPhi/.{opts}/.Options[FitMLE])===LerchPhi


LerchDistribution/: FitMLE[LerchDistribution[z1_, s1_, v1_, int_.], data1_?VectorQ, opts___?OptionQ] := Module[{e1,e2,e3,freq,cat,n,ret,i,z,s,v,ini,a,b,data}, 

   n = Length[data1];
   ini = {z1, s1, v1};

   freq = Frequency /. {opts} /. Options[FitMLE];
   If[(freq =!= True) && (freq =!= False), Message[FitMLE::type]; freq = False];

   cat = Categories /. {opts} /. Options[FitMLE];
   If[freq && VectorQ[cat] && !(Length[cat]>0 && Apply[And,Map[IntegerQ,cat]] && Apply[And,Map[NonNegative,cat]]), Message[FitMLE::cat]; Return[$Failed]];

   If[freq && Length[cat]>0 && Length[cat]<n, Message[FitMLE::len]; Return[$Failed]];

   If[freq, data = data1/Apply[Plus,data1],data=data1];

   If[freq, e1 = Apply[Plus, data cat]; e2 = (Apply[Plus, data Log[#+cat]])&; e3 = (Apply[Plus, data/(#+cat)])&, e1 = Apply[Plus, data]/n; e2 = (Apply[Plus, Log[#+data]]/n)&; e3 = (Apply[Plus, 1/(#+data)]/n)&];

   a = Min[int];  b = Max[int];
   If[freq && (a>Min[cat] || b<Max[cat]), Message[FitMLE::catf]; Return[$Failed]];
   
   rhs1[z_,s_,v_] := z/(z^a Global`LerchPhiCNCT[z,s,v+a,opts]) (z^a (a/z Global`LerchPhiCNCT[z,s,v+a,opts]+Global`LerchPhiCNCT[z,s-1,v+a+1,opts]-v Global`LerchPhiCNCT[z,s,v+a+1,opts])) /; b==Infinity;

   rhs1[z_,s_,v_] := z/(z^a Global`LerchPhiCNCT[z,s,v+a,opts]-z^(b+1) Global`LerchPhiCNCT[z,s,v+b+1,opts]) (z^a (a/z Global`LerchPhiCNCT[z,s,v+a,opts]+Global`LerchPhiCNCT[z,s-1,v+a+1,opts]-v LerchPhi[z,s,v+a+1,opts]) - z^(b+1) ((b+1)/z Global`LerchPhiCNCT[z,s,v+b+1,opts]+Global`LerchPhiCNCT[z,s-1,v+b+2,opts]-v Global`LerchPhiCNCT[z,s,v+b+2,opts])) /; b!=Infinity;

   rhs2[z_,s_,v_] := Sum[z^i Log[v+a+i]/(v+a+i)^s,{i,0,Infinity}]/ Global`LerchPhiCNCT[z,s,v+a,opts] /; b==Infinity;

   rhs2[z_,s_,v_] := (z^a Sum[z^i Log[v+a+i]/(v+a+i)^s,{i,0,Infinity}]-z^(b+1) Sum[z^i Log[v+b+i+1]/(v+b+i+1)^s,{i,0,Infinity}])/(z^a Global`LerchPhiCNCT[z,s,v+a,opts]-z^(b+1) Global`LerchPhiCNCT[z,s,v+b+1,opts]) /; b!=Infinity;

   rhs3[z_,s_,v_] := Global`LerchPhiCNCT[z,s+1,v+a,opts]/Global`LerchPhiCNCT[z,s,v+a,opts] /; b==Infinity;

   rhs3[z_,s_,v_] := (z^a Global`LerchPhiCNCT[z,s+1,v+a,opts]-z^(b+1)/s Global`LerchPhiCNCT[z,s+1,v+b+1,opts])/(z^a Global`LerchPhiCNCT[z,s,v+a,opts]-z^(b+1) Global`LerchPhiCNCT[z,s,v+b+1,opts]) /; b!=Infinity;

(*   rhs1 = z (LerchPhi[z,s-1,v+1] - v LerchPhi[z,s,v+1])/LerchPhi[z,s,v];
   rhs2 = Sum[z^i Log[v+i] / (v+i)^s, {i,0,Infinity}]/LerchPhi[z,s,v];
   rhs3 = LerchPhi[z,s+1,v]/LerchPhi[z,s,v]; *)

   With[{optf = FilterOptions[FindRoot, Sequence@@Flatten[{opts}]]}, ret = {z,s,v}/.FindRoot[{e1 == rhs1[z,s,v], e2[v] == rhs2[z,s,v], e3[v] == rhs3[z,s,v]}, {z,{ini[[1]],ini[[1]] 0.9}}, {s,{ini[[2]],ini[[2]] 0.9}}, {v,{ini[[3]],ini[[3]] 0.9}}, optf]];

   Return[ret];] /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z1,s1,v1,int],opts] && (MethodPhi/.{opts}/.Options[FitMLE])===Global`LerchPhiCNCT

FitMLE::type = "Warning: option Frequency is not specified, using False."
FitMLE::cat = "Option Categories is not a vector of nonnegative integers."
FitMLE::len = "Lengths of data and Categories are incompatible."
FitMLE::catf = "List of categories is incompatible with the domain of Lerch distribution."




LerchDistribution/: FitPearsonChiSquared[LerchDistribution[z1_, s1_, v1_, int_.], data_?VectorQ, 
    acuts_?VectorQ, opts___?OptionQ] := 
  Module[{sumrng, freqs, i, n, ini, res, z, s, v, a, b},
    
    a = Min[int]; b = Max[int];
    If[a > Min[acuts] || b < Max[acuts], Message[FitPearsonChiSquared::catf]; 
      Return[$Failed]];
    
    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]];
    
    n = Length[data];
    ini = {z1, s1, v1};
    
    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]];
    
    freq = Frequency /. {opts} /. Options[FitPearsonChiSquared];
   If[(freq =!= True) && (freq =!= False), Message[FitPearsonChiSquared::type]; freq = False];
   cat = Categories /. {opts} /. Options[FitPearsonChiSquared];
   If[freq && VectorQ[cat] && !(Length[cat]>0 && Apply[And,Map[IntegerQ,cat]] && Apply[And,Map[NonNegative,cat]]), Message[FitPearsonChiSquared::cat]; Return[$Failed]];
   If[freq && Length[cat]>0 && Length[cat]<n, Message[FitPearsonChiSquared::len]; Return[$Failed]];
   If[freq && (a>Min[cat] || b<Max[cat]), Message[FitPearsonChiSquared::catc]; Return[$Failed]];
        
    sumrng = Map[{#[[1]], #[[2]] - 1} &, pairs[acuts, 2]];
    If[freq, freqs=putfreq[data,cat,acuts], freqs = RangeCounts[data, Take[acuts, {2, Length[acuts] - 1}]]];

   If[freq, n=Apply[Plus,data]];
       
    model[z_, s_, v_] := 
      Module[{theor}, 
        theor = Map[
            Sum[PDF[LerchDistribution[z, s, v, int], i, opts], 
                Evaluate[Join[{i}, sumrng[[#]]]]] &, 
            Range[1, Dimensions[sumrng][[1]]]]; 
        n Apply[Plus, (theor - freqs/n)^2/theor]];
    
    With[{optfm = FilterOptions[FindMinimum, Sequence @@ Flatten[{opts}]]}, 
      res = FindMinimum[
          model[z, s, 
            v], {z, {ini[[1]], ini[[1]] 0.9}}, {s, {ini[[2]], 
              ini[[2]] 0.9}}, {v, {ini[[3]], ini[[3]] 0.9}}, optfm]];
    
    Return[{z, s, v} /. res[[2]]];] /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z1,s1,v1,int],opts] && ((MethodPhi/.{opts}/.Options[FitPearsonChiSquared])===Global`LerchPhiCNCT || (MethodPhi/.{opts}/.Options[FitPearsonChiSquared])===LerchPhi) && 0<z1<1 && Min[int]+v1>0 && v1!=0


LerchDistribution/: FitPearsonChiSquared[LerchDistribution[z1_, s1_, 0, int_.], data_?VectorQ, 
    acuts_?VectorQ, opts___?OptionQ] := 
  Module[{sumrng, freqs, i, n, ini, res, z, s, v, a, b},
    
    a = Min[int]; b = Max[int];
    If[a > Min[acuts] || b < Max[acuts], Message[FitPearsonChiSquared::catf]; 
      Return[$Failed]];
    
    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]];
    
    n = Length[data];
    ini = {z1, s1};
    
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]];
    
    freq = Frequency /. {opts} /. Options[FitPearsonChiSquared];
   If[(freq =!= True) && (freq =!= False), Message[FitPearsonChiSquared::type]; freq = False];
   cat = Categories /. {opts} /. Options[FitPearsonChiSquared];
   If[freq && VectorQ[cat] && !(Length[cat]>0 && Apply[And,Map[IntegerQ,cat]] && Apply[And,Map[NonNegative,cat]]), Message[FitPearsonChiSquared::cat]; Return[$Failed]];
   If[freq && Length[cat]>0 && Length[cat]<n, Message[FitPearsonChiSquared::len]; Return[$Failed]];
   If[freq && (a>Min[cat] || b<Max[cat]), Message[FitPearsonChiSquared::catc]; Return[$Failed]];
        
    sumrng = Map[{#[[1]], #[[2]] - 1} &, pairs[acuts, 2]];
    If[freq, freqs=putfreq[data,cat,acuts], freqs = RangeCounts[data, Take[acuts, {2, Length[acuts] - 1}]]];

   If[freq, n=Apply[Plus,data]];
    
    model[z_, s_, v_] := 
      Module[{theor}, 
        theor = Map[
            Sum[PDF[LerchDistribution[z, s, v, int], i, opts], 
                Evaluate[Join[{i}, sumrng[[#]]]]] &, 
            Range[1, Dimensions[sumrng][[1]]]]; 
        n Apply[Plus, (theor - freqs/n)^2/theor]];
    
    With[{optfm = FilterOptions[FindMinimum, Sequence @@ Flatten[{opts}]]}, 
      res = FindMinimum[
          model[z, s, 
            0], {z, {z1, z1 0.9}}, {s, {s1, 
              s1 0.9}}, optfm]];
    
    Return[{z, s, 0} /. res[[2]]];] /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z1,s1,0,int],opts] && (MethodPhi/.{opts}/.Options[FitPearsonChiSquared])===LerchPhi && 0<z1<1


LerchDistribution/: FitPearsonChiSquared[LerchDistribution[1, s1_, 0, int_.], data_?VectorQ, 
    acuts_?VectorQ, opts___?OptionQ] := 
  Module[{sumrng, freqs, i, n, ini, res, z, s, v, a, b},
    
    a = Min[int]; b = Max[int];
    If[a > Min[acuts] || b < Max[acuts], Message[FitPearsonChiSquared::catf]; 
      Return[$Failed]];
    
    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]];
    
    n = Length[data];
    ini = {s1};
    
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]];
    
    freq = Frequency /. {opts} /. Options[FitPearsonChiSquared];
   If[(freq =!= True) && (freq =!= False), Message[FitPearsonChiSquared::type]; freq = False];
   cat = Categories /. {opts} /. Options[FitPearsonChiSquared];
   If[freq && VectorQ[cat] && !(Length[cat]>0 && Apply[And,Map[IntegerQ,cat]] && Apply[And,Map[NonNegative,cat]]), Message[FitPearsonChiSquared::cat]; Return[$Failed]];
   If[freq && Length[cat]>0 && Length[cat]<n, Message[FitPearsonChiSquared::len]; Return[$Failed]];
   If[freq && (a>Min[cat] || b<Max[cat]), Message[FitPearsonChiSquared::catc]; Return[$Failed]];
        
    sumrng = Map[{#[[1]], #[[2]] - 1} &, pairs[acuts, 2]];
    If[freq, freqs=putfreq[data,cat,acuts], freqs = RangeCounts[data, Take[acuts, {2, Length[acuts] - 1}]]];

   If[freq, n=Apply[Plus,data]];
    
    model[z_, s_, v_] := 
      Module[{theor}, 
        theor = Map[
            Sum[PDF[LerchDistribution[z, s, v, int], i, opts], 
                Evaluate[Join[{i}, sumrng[[#]]]]] &, 
            Range[1, Dimensions[sumrng][[1]]]]; 
        n Apply[Plus, (theor - freqs/n)^2/theor]];
    
    With[{optfm = FilterOptions[FindMinimum, Sequence @@ Flatten[{opts}]]}, 
      res = FindMinimum[
          model[1, s, 
            0], {s, {s1, 
              s1 0.9}}, optfm]];
    
    Return[{1, s, 0} /. res[[2]]];] /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[1,s1,0,int],opts] && (MethodPhi/.{opts}/.Options[FitPearsonChiSquared])===LerchPhi


LerchDistribution/: FitPearsonChiSquared[LerchDistribution[1, s1_, v1_, int_.], data_?VectorQ, 
    acuts_?VectorQ, opts___?OptionQ] := 
  Module[{sumrng, freqs, i, n, ini, res, z, s, v, a, b},
    
    a = Min[int]; b = Max[int];
    If[a > Min[acuts] || b < Max[acuts], Message[FitPearsonChiSquared::catf]; 
      Return[$Failed]];
    
    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]];
    
    n = Length[data];
    ini = {s1, v1};
    
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]];
    
    freq = Frequency /. {opts} /. Options[FitPearsonChiSquared];
   If[(freq =!= True) && (freq =!= False), Message[FitPearsonChiSquared::type]; freq = False];
   cat = Categories /. {opts} /. Options[FitPearsonChiSquared];
   If[freq && VectorQ[cat] && !(Length[cat]>0 && Apply[And,Map[IntegerQ,cat]] && Apply[And,Map[NonNegative,cat]]), Message[FitPearsonChiSquared::cat]; Return[$Failed]];
   If[freq && Length[cat]>0 && Length[cat]<n, Message[FitPearsonChiSquared::len]; Return[$Failed]];
   If[freq && (a>Min[cat] || b<Max[cat]), Message[FitPearsonChiSquared::catc]; Return[$Failed]];
        
    sumrng = Map[{#[[1]], #[[2]] - 1} &, pairs[acuts, 2]];
    If[freq, freqs=putfreq[data,cat,acuts], freqs = RangeCounts[data, Take[acuts, {2, Length[acuts] - 1}]]];

   If[freq, n=Apply[Plus,data]];
    
    model[z_, s_, v_] := 
      Module[{theor}, 
        theor = Map[
            Sum[PDF[LerchDistribution[z, s, v, int], i, opts], 
                Evaluate[Join[{i}, sumrng[[#]]]]] &, 
            Range[1, Dimensions[sumrng][[1]]]]; 
        n Apply[Plus, (theor - freqs/n)^2/theor]];
    
    With[{optfm = FilterOptions[FindMinimum, Sequence @@ Flatten[{opts}]]}, 
      res = FindMinimum[
          model[1, s, 
            v], {s, {s1, s1 0.9}}, {v, {v1, 
              v1 0.9}}, optfm]];
    
    Return[{1, s, v} /. res[[2]]];] /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[1,s1,v1,int],opts] && (MethodPhi/.{opts}/.Options[FitPearsonChiSquared])===LerchPhi && Min[int]+v1>0

FitPearsonChiSquared::type = "Warning: option Frequency is not specified, using False."
FitPearsonChiSquared::cat = "Option Categories is not a vector of nonnegative integers."
FitPearsonChiSquared::len = "Lengths of data and Categories are incompatible."    
FitPearsonChiSquared::catf = "List of cutoffs is incompatible with the domain of Lerch distribution."
FitPearsonChiSquared::catc = "List of categories is incompatible with the domain of Lerch distribution."



AdjustedGeneralizedPoissonDistribution/: PDF[AdjustedGeneralizedPoissonDistribution[theta_, alpha_, phi_], x_, 
    opts___?OptionQ] := 
  If[x == 0, 
    phi + (1 - phi) Exp[-theta], (1 - phi) (1 + alpha x)^(x - 1)/
        Factorial[x] (theta Exp[-alpha theta])^x/Exp[theta]]
        
        
GeneralizedPoissonDistribution/: PDF[GeneralizedPoissonDistribution[theta_, lambda_], x_, 
    opts___?OptionQ] := 
  If[lambda < 0 && x > IntegerPart[-theta/lambda], 0, theta (theta + lambda x)^(x - 1) Exp[-theta - lambda x] / Factorial[x]]



PearsonChiSquaredStatistic[dist_, 
    data_?VectorQ, acuts_?VectorQ, opts___?OptionQ] := 
  Module[{sumrng, freqs, i, n, a, b},
    
    a = Min[int]; b = Max[int];
    If[a > Min[acuts] || b < Max[acuts], Message[PearsonChiSquaredStatistic::catf]; 
      Return[$Failed]];
    
    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]];
    
    n = Length[data];
    
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]];
    
    freq = Frequency /. {opts} /. Options[FitPearsonChiSquared];
   If[(freq =!= True) && (freq =!= False), Message[FitPearsonChiSquared::type]; freq = False];
   cat = Categories /. {opts} /. Options[FitPearsonChiSquared];
   If[freq && VectorQ[cat] && !(Length[cat]>0 && Apply[And,Map[IntegerQ,cat]] && Apply[And,Map[NonNegative,cat]]), Message[FitPearsonChiSquared::cat]; Return[$Failed]];
   If[freq && Length[cat]>0 && Length[cat]<n, Message[FitPearsonChiSquared::len]; Return[$Failed]];
   If[freq && (a>Min[cat] || b<Max[cat]), Message[FitPearsonChiSquared::catc]; Return[$Failed]];
        
    sumrng = Map[{#[[1]], #[[2]] - 1} &, pairs[acuts, 2]];
    If[freq, freqs=putfreq[data,cat,acuts], freqs = RangeCounts[data, Take[acuts, {2, Length[acuts] - 1}]]];

   If[freq, n=Apply[Plus,data]];
    
    model[z_, s_, v_] := 
      Module[{theor}, 
        theor = Map[
            Sum[PDF[dist, i, opts], 
                Evaluate[Join[{i}, sumrng[[#]]]]] &, 
            Range[1, Dimensions[sumrng][[1]]]]; 
        n Apply[Plus, (theor - freqs/n)^2/theor]];
    
    Return[model[z1, s1, v1]]] 
    
    (*/; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z1,s1,v1,int],opts] && ((MethodPhi/.{opts}/.Options[PearsonChiSquaredStatistic])===Global`LerchPhiCNCT || (MethodPhi/.{opts}/.Options[PearsonChiSquaredStatistic])===LerchPhi) *)

PearsonChiSquaredStatistic::catf = "List of cutoffs is incompatible with the domaiin of Lerch distribution."



LerchDistribution/: MMVarianceCovariance[LerchDistribution[z_,s_,v_,int_.], opts___?OptionQ] := Module[{h11,h12,h13,h21,h22,h23,h31,h32,h33,vm,h,m1,m2,m3,m4,m5,m6, l0,l5,l4,l1,l2,l11,l7,l6,l9,l8, s0,s1,s2,s3, i, vt, var},  

m1 = UncorrectedMoment[LerchDistribution[z,s,v,int],1,opts];
m2 = UncorrectedMoment[LerchDistribution[z,s,v,int],2,opts];
m3 = UncorrectedMoment[LerchDistribution[z,s,v,int],3,opts];
m4 = UncorrectedMoment[LerchDistribution[z,s,v,int],4,opts];
m5 = UncorrectedMoment[LerchDistribution[z,s,v,int],5,opts];
m6 = UncorrectedMoment[LerchDistribution[z,s,v,int],6,opts];

vm = {m2-m1^2, m4-m2^2, m6-m3^2, m3-m1 m2, m4-m1 m3, m5-m2 m3};

l0 = LerchPhi[z,s,v];
l5 = LerchPhi[z,s-2,v+1];
l4 = LerchPhi[z,s-1,v+1];
l1 = LerchPhi[z,s-1,v];
l2 = LerchPhi[z,s,v+1];
l11 = LerchPhi[z,s+1,v];
l7 = LerchPhi[z,s-3,v+1];
l6 = LerchPhi[z,s-2,v];
l9 = LerchPhi[z,s-4,v+1];
l8 = LerchPhi[z,s-3,v];
s0 = NSum[Log[v+i] z^i/(v+i)^s,{i,0,Infinity}];
s1 = NSum[Log[v+i] z^i/(v+i)^(s-1),{i,0,Infinity}];
s2 = NSum[Log[v+i] z^i/(v+i)^(s-2),{i,0,Infinity}];
s3 = NSum[Log[v+i] z^i/(v+i)^(s-3),{i,0,Infinity}];

h11 = (l5 l0-l4 l1+v (l2 l1-l4 l0))/l0^2;
h12 = (l1 s0-l0 s1)/l0^2;
h13 = ((1-s) l0^2 + s l1 l11)/l0^2 - 1;
h21 = -2 v h11 + (l7 l0-l4 l6+v (l2 l6-l5 l0))/l0^2;
h22 = -2 v h12 + (l6 s0-l0 s2)/l0^2;
h23 = -2 v h13 - 2 l1/l0 + ((2-s) l1 l0+s l11 l6)/l0^2;
h31 = -3 v^2 h11-3 v h21 + (l9 l0-l4 l8+v (l2 l8-l7 l0))/l0^2;
h32 = -3 v^2 h12-3 v h22 + (l8 s0-l0 s3)/l0^2;
h33 = -3 v^2 h13-3 v h23-3 l6/l0 + ((3-s) l6 l0+s l11 l8)/l0^2;

h = {
{h11^2, h12^2, h13^2, 2 h11 h12, 2 h11 h13, 2 h12 h13}, 
{h21^2, h22^2, h23^2, 2 h21 h22, 2 h21 h23, 2 h22 h23},
{h31^2, h32^2, h33^2, 2 h31 h32, 2 h31 h33, 2 h32 h33},
{h11 h21, h12 h22, h13 h23, h11 h22+h12 h21, h11 h23+h13 h21, h12 h23+h13 h22},
{h11 h31, h12 h32, h13 h33, h11 h32+h12 h31, h11 h33+h13 h31, h12 h33+h13 h32},
{h21 h31, h22 h32, h23 h33, h21 h32+h22 h31, h21 h33+h23 h31, h22 h33+h23 h32}};

vt = Inverse[h].vm;

var = {{vt[[1]],vt[[4]],vt[[5]]}, {vt[[4]],vt[[2]],vt[[6]]}, {vt[[5]],vt[[6]],vt[[3]]}};

Return[var];] /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z,s,v,int],opts] && (MethodPhi/.{opts}/.Options[MMVarianceCovariance])===LerchPhi


LerchDistribution/: MMVarianceCovariance[LerchDistribution[z_,s_,v_,int_.], opts___?OptionQ] := Module[{h11,h12,h13,h21,h22,h23,h31,h32,h33,vm,h,m1,m2,m3,m4,m5,m6, l0,l5,l4,l1,l2,l11,l7,l6,l9,l8, s0,s1,s2,s3, i, vt, var},  

m1 = UncorrectedMoment[LerchDistribution[z,s,v,int],1,opts];
m2 = UncorrectedMoment[LerchDistribution[z,s,v,int],2,opts];
m3 = UncorrectedMoment[LerchDistribution[z,s,v,int],3,opts];
m4 = UncorrectedMoment[LerchDistribution[z,s,v,int],4,opts];
m5 = UncorrectedMoment[LerchDistribution[z,s,v,int],5,opts];
m6 = UncorrectedMoment[LerchDistribution[z,s,v,int],6,opts];

vm = {m2-m1^2, m4-m2^2, m6-m3^2, m3-m1 m2, m4-m1 m3, m5-m2 m3};

l0 = Global`LerchPhiCNCT[z,s,v,opts];
l5 = Global`LerchPhiCNCT[z,s-2,v+1,opts];
l4 = Global`LerchPhiCNCT[z,s-1,v+1,opts];
l1 = Global`LerchPhiCNCT[z,s-1,v,opts];
l2 = Global`LerchPhiCNCT[z,s,v+1,opts];
l11 = Global`LerchPhiCNCT[z,s+1,v,opts];
l7 = Global`LerchPhiCNCT[z,s-3,v+1,opts];
l6 = Global`LerchPhiCNCT[z,s-2,v,opts];
l9 = Global`LerchPhiCNCT[z,s-4,v+1,opts];
l8 = Global`LerchPhiCNCT[z,s-3,v,opts];
s0 = NSum[Log[v+i] z^i/(v+i)^s,{i,0,Infinity}];
s1 = NSum[Log[v+i] z^i/(v+i)^(s-1),{i,0,Infinity}];
s2 = NSum[Log[v+i] z^i/(v+i)^(s-2),{i,0,Infinity}];
s3 = NSum[Log[v+i] z^i/(v+i)^(s-3),{i,0,Infinity}];

h11 = (l5 l0-l4 l1+v (l2 l1-l4 l0))/l0^2;
h12 = (l1 s0-l0 s1)/l0^2;
h13 = ((1-s) l0^2 + s l1 l11)/l0^2 - 1;
h21 = -2 v h11 + (l7 l0-l4 l6+v (l2 l6-l5 l0))/l0^2;
h22 = -2 v h12 + (l6 s0-l0 s2)/l0^2;
h23 = -2 v h13 - 2 l1/l0 + ((2-s) l1 l0+s l11 l6)/l0^2;
h31 = -3 v^2 h11-3 v h21 + (l9 l0-l4 l8+v (l2 l8-l7 l0))/l0^2;
h32 = -3 v^2 h12-3 v h22 + (l8 s0-l0 s3)/l0^2;
h33 = -3 v^2 h13-3 v h23-3 l6/l0 + ((3-s) l6 l0+s l11 l8)/l0^2;

h = {
{h11^2, h12^2, h13^2, 2 h11 h12, 2 h11 h13, 2 h12 h13}, 
{h21^2, h22^2, h23^2, 2 h21 h22, 2 h21 h23, 2 h22 h23},
{h31^2, h32^2, h33^2, 2 h31 h32, 2 h31 h33, 2 h32 h33},
{h11 h21, h12 h22, h13 h23, h11 h22+h12 h21, h11 h23+h13 h21, h12 h23+h13 h22},
{h11 h31, h12 h32, h13 h33, h11 h32+h12 h31, h11 h33+h13 h31, h12 h33+h13 h32},
{h21 h31, h22 h32, h23 h33, h21 h32+h22 h31, h21 h33+h23 h31, h22 h33+h23 h32}};

vt = Inverse[h].vm;

var = {{vt[[1]],vt[[4]],vt[[5]]}, {vt[[4]],vt[[2]],vt[[6]]}, {vt[[5]],vt[[6]],vt[[3]]}};

Return[var];] /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z,s,v,int],opts] && (MethodPhi/.{opts}/.Options[MMVarianceCovariance])===Global`LerchPhiCNCT




LerchDistribution/: MLEVarianceCovariance[LerchDistribution[z_, s_, v_, int_.], opts___?OptionQ] := Module[{i11,i12,i13,i22,i23,i33,l0,l1,l2,l3,l4,l5,l6,l7,l8,l9, s0,s1,s2,s3, i,n},  

l0 = LerchPhi[z,s,v];
l1 = LerchPhi[z,s-2,v+2];
l2 = LerchPhi[z,s-1,v+2];
l3 = LerchPhi[z,s,v+2];
l4 = LerchPhi[z,s-1,v];
l5 = LerchPhi[z,s,v+1];
l6 = LerchPhi[z,s-1,v+1];
l7 = LerchPhi[z,s+1,v+1];
l8 = LerchPhi[z,s+1,v];
l9 = LerchPhi[z,s+2,v];

s0 = NSum[Log[v+n] z^n/(v+n)^s,{n,0,Infinity}];
s1 = NSum[Log[v+n+1] (n+1) z^n/(v+n+1)^s,{n,0,Infinity}];
s2 = NSum[(Log[v+n])^2 z^n/(v+n)^s,{n,0,Infinity}];
s3 = NSum[(s Log[v+n]-1) z^n/(v+n)^(s+1),{n,0,Infinity}];

i11 = (l0 (l1-(2 v+1) l2+v (v+1) l3) - (l4-v l5)^2)/l0^2 + 1/z^2 (l4/l0-v);
i12 = ((l6-v l5) s0 - l0 s1)/l0^2;
i13 = (l0 (v l7-l5)+l8 (l6-v l5))/l0^2 s;
i22 = (l0 s2 - s0^2)/l0^2;
i23 = (l0 s3 - s l8 s0)/l0^2 + l8/l0;
i33 = (l0 l9-l8^2)/l0^2 s^2;

i = {{i11,i12,i13},{i12,i22,i23},{i13,i23,i33}};

Return[Inverse[i]];]  /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z,s,v,int],opts]  && (MethodPhi/.{opts}/.Options[MLEVarianceCovariance])===LerchPhi


LerchDistribution/: MLEVarianceCovariance[LerchDistribution[z_, s_, v_, int_.], opts___?OptionQ] := Module[{i11,i12,i13,i22,i23,i33,l0,l1,l2,l3,l4,l5,l6,l7,l8,l9, s0,s1,s2,s3, i,n},  

l0 = Global`LerchPhiCNCT[z,s,v,opts];
l1 = Global`LerchPhiCNCT[z,s-2,v+2,opts];
l2 = Global`LerchPhiCNCT[z,s-1,v+2,opts];
l3 = Global`LerchPhiCNCT[z,s,v+2,opts];
l4 = Global`LerchPhiCNCT[z,s-1,v,opts];
l5 = Global`LerchPhiCNCT[z,s,v+1,opts];
l6 = Global`LerchPhiCNCT[z,s-1,v+1,opts];
l7 = Global`LerchPhiCNCT[z,s+1,v+1,opts];
l8 = Global`LerchPhiCNCT[z,s+1,v,opts];
l9 = Global`LerchPhiCNCT[z,s+2,v,opts];

s0 = NSum[Log[v+n] z^n/(v+n)^s,{n,0,Infinity}];
s1 = NSum[Log[v+n+1] (n+1) z^n/(v+n+1)^s,{n,0,Infinity}];
s2 = NSum[(Log[v+n])^2 z^n/(v+n)^s,{n,0,Infinity}];
s3 = NSum[(s Log[v+n]-1) z^n/(v+n)^(s+1),{n,0,Infinity}];

i11 = (l0 (l1-(2 v+1) l2+v (v+1) l3) - (l4-v l5)^2)/l0^2 + 1/z^2 (l4/l0-v);
i12 = ((l6-v l5) s0 - l0 s1)/l0^2;
i13 = (l0 (v l7-l5)+l8 (l6-v l5))/l0^2 s;
i22 = (l0 s2 - s0^2)/l0^2;
i23 = (l0 s3 - s l8 s0)/l0^2 + l8/l0;
i33 = (l0 l9-l8^2)/l0^2 s^2;

i = {{i11,i12,i13},{i12,i22,i23},{i13,i23,i33}};

Return[Inverse[i]];]  /; MatchQ[int, HoldPattern[Interval[{_?(Or[Positive[#], # == 0] &), _?(Or[Positive[#], # == 0] &)}]]] && ParameterQ[LerchDistribution[z,s,v,int],opts]  && (MethodPhi/.{opts}/.Options[MLEVarianceCovariance])===Global`LerchPhiCNCT



Protect @@ prt

End[]
EndPackage[]