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

Of the Mathematica package SDistribution.m

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


This software (file SDistribution.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 authors 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: 0.99 *)

(*:Name: SDistribution *)

(*:Context: SDistribution`SDistribution` *)

(*:Title: Calculation and fitting of the S distribution *)

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

(*:History:
   Created: August, 2001 by Sergej V. Aksenov.   
   Last modified: May, 2002 by Sergej V. Aksenov.*)

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

(*:Reference: 
Usage messages, documentation notebook SDistribution.nb, article "Mathematica and C programs for minimum distance estimation of the S distribution and for calculation of goodness-of-fit by bootstrap", by Sergej V. Aksenov and Michael A. Savageau, submitted to Statistics and Computing (2002). *)

(*:Summary: 
This package provides Mathematica functions that define the S distribution and calculate its properties and estimate its parameters. *)

(*:Keywords:  Continuous distribution, S distribution, parameter estimation, minimum distance, Lerch's transcendent, optimization *)

(*:Limitations: 
*  Moments are not implemented (this work is in progress; I plan to add Mean, Variance, Standarddeviation, Skewness, and Kurtosis calculated by Voit's ode method..

*  Diagnostic printout from external libraries CVODE, DFO, CFSQo is not used; should be easy to write it to files in the working directory..  a To Do item..
*)

(*:Note:



***  TO DO:
use functional substitution instead of transformation for v<0 in quantile
See line 604 below
***



*  The following standard functions have been modified to work with the S distribution: ParameterQ, Domain, DomainQ, Quantile, Random, RandomArray, PDF, CDF. 

*  When specifying parameters for the S distribution, care must be taken when calculating quantiles.  If g and h are such that g=(h*k+1)/(k+1) for integer k, there will be loss of precision when calculating k=(g-1)/(h-g) in the routine that calculates quantiles.  This routine needs the correct k to determine which method to use.  It appears as though subtracting in the denominator and numerator of this formula loses precision; therefore a temporary hack is that k=(g*10-1*10)/(h*10-g*10).  Alternatively, one might want to give values of g and h in extended precision to avoid these cancellation errors. 
 
*  SetOptions[DomainQ, LerchPhiRelativeAccuracy->user value, MethodPhi->...] will not take effect; instead apply SetOptions to Domain, which is called by DomainQ.  Also, whatever LerchPhiRelativeAccuracy or MethodPhi is set for Domain/DomainQ, it will be overriden by a corresponding setting of a function that calls Domain/DomainQ (e.g., PDF, CDF).

*  Options of standard functions are reinitialized here; they don't seem to have options in a standard setting  so this should not be a problem, but...  may be careful..

*  Parameters are not always checked for proper type, or if they are numeric etc., especially those passed as options.  Therefore use all functions with care.  Installed MathLink programs may freeze or exit if supplied with arguments of a wrong type.

*  GoodnessOfFitStatistic does not check and will not complain if input SDistribution parameters are not consistent with the data (it calls CDF that does the check; CDF returns 0 if out of domain). 

*  Abort by standard Mathematica way e.g. Cmd(cntrl)-dot is implemented for optim MathLink program only; lerch and cdf don't perform loops that can become long etc.  Also, note that no yield functions are called, optim just checks for MLAbort!=0.  This works fine on Unix systems, but may not wotk on older Mac or Windows without preemptive multitasking..

*  Quantile functions are implemented also with standard Mathematica LerchPhi, however since it calculates according to "alternative" definition, wrong results obtsain for v<0, s odd (this I took care of by using transformation to positive v);  moreover, for v<0 very close to integer both Mathematica and CNCT implementations of LerchPhi lose precision and result calcuates incorrectly.  For Mathematica LerchPhi the remedy is to supply Quantle etc. functions with arguments in quadruple precision. For the binary-provided Phi nothing can be done, except to implement the C code in extended precision (which exists, just have to do it someday..)  Note however, that in optim.exe CNCT does not encounter loss of precision problems because v>0 always (g<1).
*)



BeginPackage["SDistribution`SDistribution`", "Statistics`DescriptiveStatistics`", "Statistics`Common`DistributionsCommon`", "Statistics`ContinuousDistributions`", "Statistics`DataManipulation`", "Statistics`NonlinearFit`", "Statistics`Common`RegressionCommon`", "Statistics`NormalDistribution`", "Statistics`ConfidenceIntervals`", "Utilities`FilterOptions`", "Graphics`Graphics`"]             
                          
                          
                          
Unprotect[SDistribution];



(* Function names that will be exported. *)

SDistribution::usage = 
"SDistribution[g, h, alpha, x0, f0] represents the S distribution with the shape parameters g and h, scale parameter alpha, and location parameter x0 such that F(x0)=f0."

SDistributionParameterEstimate::usage = 
"SDistributionParameterEstimate[data] returns a list of estimates of the S distribution parameters computed from data."

SDistributionOptimize::usage = 
"SDistributionOptimize[data, par] returns a list of S distribution parameters that minimizes a goodness-of-fit statistic and achieved minimum value of the statistic for data and initial estimates par."

EmpiricalDF::usage = 
"EmpiricalDF[data] returns a list {data, ecdf} of an ordered list of unique points from data and a corresponding list empirical cumulative distribution function values."

GoodnessOfFitStatistic::usage = 
"GoodnessOfFitStatistic[data, distribution] returns value of a goodness-of-fit statistic computed from data for the completely specified distribution."

SampleQuantileCI::usage = 
"SampleQuantileCI[data, q] returns a list{lower, upper} representing nonparametric confidence interval for the qth sample quantile based on fractional order statistics of data."

LerchPhiCNCT::usage =
"LerchPhi[z, s, v] gives the Lerch's transcendent Phi(z, s, v) calculated with the combined-nonlinear condensation transformation." 



(* Option names that will be exported. *)

EmpiricalDFDefinition::usage=
"Option EmpiricalDFDefinition specifies whether ties should be eliminated from data giving proportionately more weight on the unique points (EliminateTies) or tied points should be kept (KeepTies)."

EliminateTies::usage=
"Value for option EmpiricalDFDefinition that causes to eliminate tied observations."

KeepTies::usage=
"Value for option EmpiricalDFDefinition that causes to keep tied observations."

EmpiricalDFDefinitionConstant::usage=
"Option EmpiricalDFDefinitionConstant specifies the constant in the formula for empirical d.f. if the ties are kept."

SampleQuantile::usage =
"Option SampleQuantile specifies the sample quantile that will be used to estimate S distribution parameter x0."

EstimateMethod::usage =
"Option EstimateMethod specifies method to use in obtaining estimates of S distribution parameters."

NumberPointsInBins::usage =
"Option NumberPointsInBins specifies the number of data points per bin from which the number of bins is calculated."

AlgebraicRegression::usage =
"Method of estimating parameters of the S distribution."

QuantileEquation::usage =
"Method of estimating parameters of the S distribution."

ApproximateMLE::usage =
"Method of estimating parameters of the S distribution."

LerchPhiRelativeAccuracy::usage =
"Option LerchPhiRelativeAccuracy specifies digits of relative accuracy of Lerch's transcendent."

ZeroinIntervalUncertainty::usage =
"Option ZeroinIntervalUncertainty specifies the length of the interval of uncertainty of the final result in the zero finding routine."

CVODEAbsoluteTolerance::usage =
"Option CVODEAbsoluteTolerance specifies digits of the absolute accuracy of final result in the CVODE library."

CVODERelativeTolerance::usage =
"Option CVODERelativeTolerance specifies digits of the relative accuracy of final result in the CVODE library."

CVODEIntegrationMethod::usage =
"Option CVODEIntegrationMethod specifies the type of linear multistep method that CVODE library uses to integrate differential equations (BDF or ADAMS)."

BDF::usage =
"BDF is value for option CVODEIntegrationMethod that specifies BDF formula as the linear multistep method."

ADAMS::usage =
"ADAMS is value for option CVODEIntegrationMethod that specifies Adams-Moulton formula as the linear multistep method."

CVODEIterationType::usage =
"Option CVODEIterationType specifies the type of iteration used to solve the nonlinear system that arises during each internal time step in the CVODE library (NEWTON or FUNCTIONAL)."

NEWTON::usage =
"NEWTON is value for option CVODEIterationType that specifies the functional iteration for solution of the nonlinear system at each time step."

FUNCTIONAL::usage =
"FUNCTIONAL is value for option CVODEIterationType that specifies Newton iteration for solution of the nonlinear system at each time step."

CVODEJacobian::usage =
"Option CVODEJacobian specifies how the Jacobian of the system of differential equations is computed in the CVODE library (ANALYTICAL or NUMERICAL)."

ANALYTICAL::usage =
"ANALYTICAL is value for option CVODEJacobian that specifies the user-provided jacobian in a linear solver module."

NUMERICAL::usage =
"NUMERICAL is value for option CVODEJacobian that specifies the difference quotient approximation for the jacobian in a linear solver module."

CVODEMaxSteps::usage =
"Option CVODEMaxSteps specifies the maximum number of internal steps to be taken by the CVODE solver in its attempt to reach tout."

DFOTrustLowerBound::usage = 
"Option DFOTrustLowerBound specifies the lower bound on the trust region (stopping criterion)."

DFOTrustUpperBound::usage = 
"Option DFOTrustUpperBound specifies the upper bound on the trust region."

DFOTrustInitial::usage = 
"Option DFOTrustInitial specifies the iinitial trust region radius."

DFOObjectiveLowerBound::usage = 
"Option DFOObjectiveLowerBound specifies the lower bound on the objective function."

DFOMaxFunctionEvaluations::usage = 
"Option DFOMaxFunctionEvaluations specifies the maximum number of function evaluations allowed by the algorithm."

DFOMaxIterations::usage =
"Option DFOMaxIterations specifies the maximum number of iterations allowed."

DFOAbsoluteFunctionNoise::usage =
"Option DFOAbsoluteFunctionNoise specifies the absolute noise of function value computations."

DFORelativeFunctionNoise::usage =
"Option DFORelativeFunctionNoise specifies the relative noise of function value computations."

DFOPivotThreshold::usage =
"Option DFOPivotThreshold specifies the threshold on the pivot size for accepting a new point into the interpolation set."

DFOPivotThresholdGeometry::usage =
"Option DFOPivotThresholdGeometry specifies the threshold on the pivot size for accepting a new point into the interpolation set, while caring about the geometry."

DFOPivotFactor::usage =
"Option DFOPivotFactor specifies the factor of improvement in the pivot, when replacing an existing interpolation point by another point, with hopes to improve the geometry."

DFOMinPointsInterpolation::usage =
"Option DFOMinPointsInterpolation specifies the minimum number of points in an interpolation."

DFOMaxBaseDistanceFactor::usage =
"Option DFOMaxBaseDistanceFactor (times DFOTrustInitial) specifies the maximum distance from the base within which the points for interpolation are selected."

DFOInterpolationEffort::usage =
"Option DFOInterpolationEffort specifies the level of effort that the algorithm should spend recomputing an interpolation set (1 -- when all else fails or the model too old, 2 -- when reduction is bad and adding new point was not successfull, 3 -- whenever reduction is not good, 4 -- every time)."

DFOQuadraticModel::usage =
"Option DFOQuadraticModel specifies the quadratic models in case the interpolation is not full (1 -- minimum Frobenius norm, 2 -- derived from Newton polynomials)."

DFOScaling::usage =
"Option DFOScaling specifies scaling (0 -- none, 1 -- based on the initial point, 2 -- based on the box constraints if present)."

(* DFOPrintoutLevel::usage =
"Option DFOPrintoutLevel specifies level of printout (-1 -- none, 0, 1, 2, 3 -- progressively more)."
DFOOutputDevice::usage =
"Option DFOOutputDevice specifies the printout device (6 prints on display)."
*)

CFSQPObjectiveAbsoluteEpsilon::usage =
"Option CFSQPObjectiveAbsoluteEpsilon specifies the requirement for the absolute difference between  values of the objective function at successive iterates (used as a stopping criterion).  Set to negative value to switch off."

CFSQPObjectiveRelativeEpsilon::usage =
"Option CFSQPObjectiveRelativeEpsilon specifies the requirement for the relative difference between  values of the objective function at successive iterates (used as a stopping criterion).  Set to negative value to switch off."

CFSQPKuhnTuckerFinalNorm::usage =
"Option CFSQPKuhnTuckerFinalNorm specifies the final norm requirement for the Kuhn-Tucker norm (used as a stopping criterion).  Set to negative value to switch off."

CFSQPMaxIterations::usage =
"Option CFSQPMaxIterations specifies maximum number of iterations for the CFSQP algorithm."

CFSQPNewtonDirectionFinalNorm::usage = 
"Option CFSQPNewtonDirectionFinalNorm specifies the final norm requirement for the Newton direction d_k^0 (eps) in the CFSQP algorithm (used as a stopping criterion).  It must be bigger than the machine precision."

CFSQPFiniteDifferencePerturbation::usage =
"Option CFSQPFiniteDifferencePerturbation specifies the perturbation size to use in approximating gradients by finite difference."

CFSQPCheckConstraintsFirst::usage =
"Option CFSQPCheckConstraintsFirst specifies whether a CFSQP algorithm at each trial point during line search should check for constraints first (True).  If False, the function (objective or constraint) that caused the previous value of t to be rejected is checked first."

CFSQPAlgorithm::usage = 
"Option CFSQPAlgorithm selects a CFSQP algorithm for optimization: FSQPAL requires objective function to decrease at each iteration, while FSQPNL requires objective function to decrease within at most four iterations."

FSQPAL::usage = 
"FSQPAL is value for option CFSQPAlgorithm that selects FSQP-AL algorithm that requires objective function to decrease at each iteration."

FSQPNL::usage = 
"FSQPNL is value for option CFSQPAlgorithm that selects FSQP-NL algorithm that requires objective function to decrease within at most four iterations."

OptimizationAlgorithm::usage = 
"Option OptimizationAlgorithm specifies algorithm to use in SDistributionOptimize (CFSQP or DFO)."

CFSQP::usage =
"CFSQP is value for option OptimizationAlgorithm."

DFO::usage =
"DFO is value for option OptimizationAlgorithm."

Statistic::usage =
"Option Statistic specifies either of six goodness-of-fit statistics: KolmogorovSmirnov, Kuiper, CramerVonMises, Watson, AndersonDarling, RMS."

KolmogorovSmirnov::usage =
"KolmogorovSmirnov is value for option Statistic."

Kuiper::usage =
"Kuiper is value for option Statistic."

CramerVonMises::usage =
"CramerVonMises is value for option Statistic."

Watson::usage =
"Watson is value for option Statistic."

AndersonDarling::usage =
"AndersonDarling is value for option Statistic."

RMS::usage =
"RMS (root mean squared error) is value for option Statistic."

MethodODE::usage =
"Option MethodODE specifies either NDSolve or CVODE as differential equation solvers."

CVODE::usage =
"CVODE is value for option MethodODE that specifies external C library CVODE as the differential equation solver in functions PDF and CDF when applied to S distributions."

MethodPhi::usage =
"Option MethodPhi specifies either LerchPhi or LerchPhiCNCT as algorithms for calculating Lerch's transcendent Phi."



(* Set default options. *)

LerchPhiDefaultOption = {LerchPhiRelativeAccuracy -> 14, MethodPhi->LerchPhiCNCT}
CVODEDefaultOptions = {MethodODE->CVODE, ZeroinIntervalUncertainty -> 0.0, CVODEAbsoluteTolerance -> 10, CVODERelativeTolerance -> 0, CVODEIntegrationMethod -> BDF, CVODEIterationType -> NEWTON, CVODEJacobian -> ANALYTICAL, CVODEMaxSteps ->1000}
SDOptimizeDefaultOptions = {Statistic -> AndersonDarling, OptimizationAlgorithm -> CFSQP, DFOTrustLowerBound -> 1.0 10^(-4), DFOTrustUpperBound -> 100., DFOTrustInitial -> 0.1, DFOObjectiveLowerBound -> -1.0 10^20, DFOMaxFunctionEvaluations -> 1500, DFOMaxIterations -> 900, DFOAbsoluteFunctionNoise -> 0.0, DFORelativeFunctionNoise -> 0.0, DFOPivotThreshold -> 0.01, DFOPivotThresholdGeometry -> 0.1, DFOPivotFactor -> 0.1, DFOMinPointsInterpolation -> 5, DFOMaxBaseDistanceFactor -> 6, DFOInterpolationEffort -> 1, DFOQuadraticModel -> 1, DFOScaling -> 0, CFSQPNewtonDirectionFinalNorm -> 1.0 10^(-5), CFSQPFiniteDifferencePerturbation -> 0.0, CFSQPAlgorithm -> FSQPNL, CFSQPMaxIterations -> 500, CFSQPCheckConstraintsFirst -> False, CFSQPObjectiveAbsoluteEpsilon -> -1.0, CFSQPObjectiveRelativeEpsilon -> -1.0, CFSQPKuhnTuckerFinalNorm -> -1.0}
(* DFOPrintoutLevel -> -1, DFOOutputDevice -> 6, *)
Options[SDistributionParameterEstimate] = {SampleQuantile -> 0.5, EstimateMethod -> AlgebraicRegression, NumberPointsInBins -> 6}

Unprotect[Quantile,Random]
Options[Quantile] = LerchPhiDefaultOption
Options[Random] = LerchPhiDefaultOption
Protect[Quantile,Random]
Options[RandomArray] = LerchPhiDefaultOption
Options[Domain] = LerchPhiDefaultOption
Options[DomainQ] = LerchPhiDefaultOption
Options[CDF] = Join[LerchPhiDefaultOption, CVODEDefaultOptions]
Options[PDF] = Join[LerchPhiDefaultOption, CVODEDefaultOptions]
Options[SDistributionOptimize] = Join[LerchPhiDefaultOption, CVODEDefaultOptions, SDOptimizeDefaultOptions]
Options[EmpiricalDF] = {EmpiricalDFDefinition -> EliminateTies, EmpiricalDFDefinitionConstant -> 0.5}
Options[GoodnessOfFitStatistic] = {Statistic -> AndersonDarling}
Options[SampleQuantileCI] = {ConfidenceLevel -> .95}
Options[LerchPhiCNCT] = LerchPhiDefaultOption


Begin["`Private`"]



(* Text of error messages. *)

SdCDFCvderr1 = "CVODE error `1`: the cvode_mem argument was NULL."
SdCDFCvderr2 = "CVODE error `1`: one of the inputs to CVODE is illegal."
SdCDFCvderr3 = "CVODE error `1`: the solver took mxstep internal steps but could not reach tout."
SdCDFCvderr4 = "CVODE error `1`: the solver could not satisfy the accuracy demanded by the user for some internal step."
SdCDFCvderr5 = "CVODE error `1`: error test failures occurred too many times during one internal time step or occurred with |h|=hmin."
SdCDFCvderr6 = "CVODE error `1`: convergence test failures occurred too many times during one internal time step or occurred with |h|=hmin."
SdCDFCvderr7 = "CVODE error `1`: the linear solver's setup routine failed in an unrecoverable manner."
SdCDFCvderr8 = "CVODE error `1`: the linear solver's solve routine failed in an unrecoverable manner."
SdCDFCvdneg = "SDCDF error: CVODE gives negative solution."

LerchDvrg1 = "LerchPhi error `1`: Lerch Phi series diverges for 1 <= |z|."
LerchIntvdef2 = "LerchPhi error `1`: Lerch Phi is not defined for integer v <= 0."
LerchAnegsr3 = "LerchPhi error `1`: C function pow() is not defined for v < 0 and s not integer in LerchPhi."
LerchIntaj4 = "LerchPhi error `1`: Integer overflow in Lerch Phi."
LerchUnderom5 = "LerchPhi error `1`: Underflow in remainder estimate omega in Lerch Phi."
LerchCnvg6 = "LerchPhi error `1`: Convergence not achieved within maximum number of iterations in Lerch Phi."

SDOptimdfo1p = "DFO termination parameter `1`: too many function evaluations."
SDOptimdfo2p = "DFO termination parameter `1`: too many iterations."
SDOptimdfo1m = "DFO termination parameter `1`: illegal value of N, NCLIN, or NCNLN."
SDOptimdfo2m = "DFO termination parameter `1`: illegal value of LDX or LDA."
SDOptimdfo3m = "DFO termination parameter `1`: real workspace is too small."
SDOptimdfo4m = "DFO termination parameter `1`: integer workspace is too small."
SDOptimdfo5m = "DFO termination parameter `1`: inconsistent bounds."
SDOptimdfo6m = "DFO termination parameter `1`: function value cannot be computed at a given initial point."
SDOptimdfo7m = "DFO termination parameter `1`: minimization over the trust region failed."
SDOptimdfo8m = "DFO termination parameter `1`: MAXNF is too small, less than 2."
SDOptimdfo9m = "DFO termination parameter `1`: cannot find a two-point model."

SDOptimcfsqp1 = "CFSQP error `1`: initial guess is infeasible for linear constraints."
SDOptimcfsqp2 = "CFSQP error `1`: initial guess is infeasible for nonlinear inequality constraints and linear constraints."
SDOptimcfsqp3 = "CFSQP error `1`: maximum number of iterations has been reached before a solution was obtained."
SDOptimcfsqp4 = "CFSQP error `1`: the line search fails to find a new iterate."
SDOptimcfsqp5 = "CFSQP error `1`: failure of the QP solver in attempting to construct d^0."
SDOptimcfsqp6 = "CFSQP error `1`: failure of the QP solver in attempting to construct d^1."
SDOptimcfsqp7 = "CFSQP error `1`: input data are not consistent."
SDOptimcfsqp8 = "CFSQP error `1`: new iterate is numerically equivalent to the previous iterate, though stopping criterion is not yet satisfied."
SDOptimcfsqp9 = "CFSQP error `1`: one of the penalty parameters exceeded bigbnd."



(* Install external binaries. *)

SDistribution::binl = 
"Installation of the MathLink binary lerch.exe failed.  Turn option LerchPhiTranscendent->BuiltIn to use Mathematica's LerchPhi."
SDistribution::binc = 
"Installation of the MathLink binary cdf.exe failed.  Turn option DifferentialEquationSolver->BuiltIn to use Mathematica's NDSolve."
SDistribution::bino = 
"Installation of the MathLink binary optim failed.  The binary must be found for the minimum distance estimation component of the package to work."

If[Install["SDistribution`lerch`"]===$Failed, Message[SDistribution::binl]]
If[Install["SDistribution`cdf`"]===$Failed, Message[SDistribution::binc]]
If[Install["SDistribution`optim`"]===$Failed, Message[SDistribution::bino]]



(* Extension of standard statistical functions to S distribution. *)

SDistribution/: ParameterQ[SDistribution[g_?NumberQ, h_?NumberQ, alpha_?NumberQ, x0_?NumberQ, f0_?NumberQ]] := And[If[N[g] < N[h], True, Message[SDistribution::ghineq, g, h]; False, True], If[N[alpha] > 0, True, Message[SDistribution::posalp, alpha]; False, True], If[0 < N[f0] && N[f0] < 1, True, Message[SDistribution::ffrac, f0]; False, True]]

SDistribution::posalp = 
"The scale parameter `1` is expected to be positive."
SDistribution::ghineq = 
"The shape parameter `1` is expected to be less than shape parameter `2`."
SDistribution::ffrac =
"Initialization cdf `1` is expected to be between 0 and 1."

SDistribution/: Domain[SDistribution[g_, h_, alpha_, x0_, f0_], opts___?OptionQ] := Module[{xasym}, If[1 <= g, Return[Interval[{-Infinity, Infinity}]], xasym = SDistributionQuantile[{g, h, alpha, x0, f0}, 0.0, 4, opts]; If[xasym===$Failed, Return[$Failed], Return[Interval[{xasym, Infinity}]]]]] /; ParameterQ[SDistribution[g, h, alpha, x0, f0]]

SDistribution/: DomainQ[SDistribution[g_, h_, alpha_, x0_, f0_], list_?VectorQ, opts___?OptionQ] := Module[{dom}, dom = Domain[SDistribution[g, h, alpha, x0, f0], FilterOptions[Domain,opts]];  If[dom===$Failed, Return[$Failed], Return[Scan[If[!IntervalMemberQ[dom, #], Return[False]]&, list] =!= False]]] /; ParameterQ[SDistribution[g, h, alpha, x0, f0]]

SDistribution/: DomainQ[SDistribution[g_, h_, alpha_, x0_, f0_], x_, opts___?OptionQ] := Module[{dom}, dom = Domain[SDistribution[g, h, alpha, x0, f0], FilterOptions[Domain,opts]]; If[dom===$Failed, Return[$Failed], Return[IntervalMemberQ[dom, x]]]] /; ParameterQ[SDistribution[g, h, alpha, x0, f0]]

SDistribution/: PDF[SDistribution[g_, h_, alpha_, x0_, f0_], x_?VectorQ, opts___?OptionQ] := Module[{F,domres,out1}, out1 = {};
Catch[
Scan[(domres = DomainQ[SDistribution[g, h, alpha, x0, f0], #, FilterOptions[DomainQ,opts]]; If[domres===$Failed, Throw[$Failed]]; If[!domres, AppendTo[out1, 0.0], F = SDistributionDF[{g, h, alpha, x0, f0}, {#}, 1, opts];  If[F=!=$Failed, AppendTo[out1, alpha (F^g - F^h)], Throw[$Failed]]])&, x]; Return[Flatten[out1]]]] /; ParameterQ[SDistribution[g, h, alpha, x0, f0]] && (MethodODE/.{opts}/.Options[PDF])===CVODE

SDistribution/: PDF[SDistribution[g_, h_, alpha_, x0_, f0_], x_?VectorQ, opts___?OptionQ] := Module[{F,domres,out1,res,xt,FF}, out1 = {};
Catch[
Scan[(domres = DomainQ[SDistribution[g, h, alpha, x0, f0], #, FilterOptions[DomainQ,opts]]; If[domres===$Failed, Throw[$Failed]]; If[!domres, AppendTo[out1, 0.0], res=NDSolve[{FF'[xt]==alpha (FF[xt]^g-FF[xt]^h),FF[x0]==f0},FF[xt],{xt,x0,#},FilterOptions[NDSolve,opts]]; F = Re[Part[FF[xt]/.res/.xt->#,1]]; If[F<0, F=0]; If[F=!=$Failed, AppendTo[out1, alpha (F^g - F^h)], Throw[$Failed]]])&, x]; Return[Flatten[out1]]]] /; ParameterQ[SDistribution[g, h, alpha, x0, f0]] && (MethodODE/.{opts}/.Options[PDF])===NDSolve

SDistribution/: PDF[SDistribution[g_, h_, alpha_, x0_, f0_], x_, opts___?OptionQ] := Module[{F,domres}, Catch[domres = DomainQ[SDistribution[g, h, alpha, x0, f0], x, FilterOptions[DomainQ,opts]]; If[domres===$Failed, Throw[$Failed]]; If[domres, F = SDistributionDF[{g, h, alpha, x0, f0}, {x}, 1, opts]; If[F=!=$Failed, Return[alpha (F[[1]]^g - F[[1]]^h)], Return[$Failed]], Return[0.]]]] /; ParameterQ[SDistribution[g, h, alpha, x0, f0]] && (MethodODE/.{opts}/.Options[PDF])===CVODE

SDistribution/: PDF[SDistribution[g_, h_, alpha_, x0_, f0_], x_, opts___?OptionQ] := Module[{F,domres,res,xt,FF}, Catch[domres = DomainQ[SDistribution[g, h, alpha, x0, f0], x, FilterOptions[DomainQ,opts]]; If[domres===$Failed, Throw[$Failed]]; If[domres, res=NDSolve[{FF'[xt]==alpha (FF[xt]^g-FF[xt]^h),FF[x0]==f0},FF[xt],{xt,x0,x},FilterOptions[NDSolve,opts]]; F = Re[Part[FF[xt]/.res/.xt->x,1]]; If[F<0, F=0]; If[F=!=$Failed, Return[alpha (F^g - F^h)], Return[$Failed]], Return[0.]]]] /; ParameterQ[SDistribution[g, h, alpha, x0, f0]] && (MethodODE/.{opts}/.Options[PDF])===NDSolve

SDistribution/: CDF[SDistribution[g_, h_, alpha_, x0_, f0_], x_?VectorQ, opts___?OptionQ] := Module[{F,domres,out1}, out1 = {};
Catch[
sc=Scan[(domres = DomainQ[SDistribution[g, h, alpha, x0, f0], #, FilterOptions[DomainQ,opts]]; If[domres===$Failed, Throw[$Failed]]; If[!domres, AppendTo[out1, 0.0], F = SDistributionDF[{g, h, alpha, x0, f0}, {#}, 2, opts]; If[F=!=$Failed, AppendTo[out1, F], Throw[$Failed]]])&, x]; Return[Flatten[out1]]]] /; ParameterQ[SDistribution[g, h, alpha, x0, f0]] && (MethodODE/.{opts}/.Options[CDF])===CVODE

SDistribution/: CDF[SDistribution[g_, h_, alpha_, x0_, f0_], x_?VectorQ, opts___?OptionQ] := Module[{F,domres,out1,res,xt,FF}, out1 = {};
Catch[
sc=Scan[(domres = DomainQ[SDistribution[g, h, alpha, x0, f0], #, FilterOptions[DomainQ,opts]]; If[domres===$Failed, Throw[$Failed]]; If[!domres, AppendTo[out1, 0.0], res=NDSolve[{FF'[xt]==alpha (FF[xt]^g-FF[xt]^h),FF[x0]==f0},FF[xt],{xt,x0,#},FilterOptions[NDSolve,opts]]; F = Re[Part[FF[xt]/.res/.xt->#,1]]; If[F<0, F=0]; If[F=!=$Failed, AppendTo[out1, F], Throw[$Failed]]])&, x]; Return[Flatten[out1]]]] /; ParameterQ[SDistribution[g, h, alpha, x0, f0]] && (MethodODE/.{opts}/.Options[CDF])===NDSolve

SDistribution/: CDF[SDistribution[g_, h_, alpha_, x0_, f0_], x_, opts___?OptionQ] := Module[{F,domres}, Catch[domres = DomainQ[SDistribution[g, h, alpha, x0, f0], x, FilterOptions[DomainQ,opts]]; If[domres===$Failed, Throw[$Failed]]; If[domres, F = SDistributionDF[{g, h, alpha, x0, f0}, {x}, 2, opts][[1]]; If[F=!=$Failed, Return[F], Return[$Failed]], Return[0.]]]] /; ParameterQ[SDistribution[g, h, alpha, x0, f0]] && (MethodODE/.{opts}/.Options[CDF])===CVODE

SDistribution/: CDF[SDistribution[g_, h_, alpha_, x0_, f0_], x_, opts___?OptionQ] := Module[{F,domres,res,xt,FF}, Catch[domres = DomainQ[SDistribution[g, h, alpha, x0, f0], x, FilterOptions[DomainQ,opts]]; If[domres===$Failed, Throw[$Failed]]; If[domres, res=NDSolve[{FF'[xt]==alpha (FF[xt]^g-FF[xt]^h),FF[x0]==f0},FF[xt],{xt,x0,x},FilterOptions[NDSolve,opts]]; F = Re[Part[FF[xt]/.res/.xt->x,1]]; If[F<0, F=0]; If[F=!=$Failed, Return[F], Return[$Failed]], Return[0.]]]] /; ParameterQ[SDistribution[g, h, alpha, x0, f0]] && (MethodODE/.{opts}/.Options[CDF])===NDSolve

SDistribution/: Quantile[SDistribution[g_, h_, alpha_, x0_, f0_], fraction_?VectorQ, opts___?OptionQ] := Module[{out}, If[!Apply[And, Map[QuantileQ[#] &, fraction]], Return[$Failed]]; out={}; Scan[Which[#==1, AppendTo[out,Infinity], (#==0&&g>=1), AppendTo[out,-Infinity], True, AppendTo[out,SDistributionQuantile[{g, h, alpha, x0, f0}, #, 1, opts]]]&, fraction]; Return[out]] /; ParameterQ[SDistribution[g, h, alpha, x0, f0]]

SDistribution/: Quantile[SDistribution[g_, h_, alpha_, x0_, f0_], fraction_, opts___?OptionQ] := Which[fraction==1, Return[Infinity], (fraction==0&&g>=1), Return[-Infinity], True, Return[SDistributionQuantile[{g, h, alpha, x0, f0}, fraction, 1, opts]]] /; QuantileQ[fraction] && ParameterQ[SDistribution[g, h, alpha, x0, f0]]

SDistribution/: Random[SDistribution[g_, h_, alpha_, x0_, f0_], opts___?OptionQ] := Module[{rr}, rr = Random[]; If[rr==1, Return[Infinity], Return[SDistributionQuantile[{g, h, alpha, x0, f0}, rr, 2, opts]]]] /; ParameterQ[SDistribution[g, h, alpha, x0, f0]]

SDistribution/: RandomArray[SDistribution[g_, h_, alpha_, x0_, f0_], dim_, opts___?OptionQ] := Module[{rr}, Table[rr = Random[]; If[rr==1, Infinity, SDistributionQuantile[{g, h, alpha, x0, f0}, rr, 3, opts]], {dim}] ] /; ParameterQ[SDistribution[g, h, alpha, x0, f0]] && ((IntegerQ[dim] && dim > 0) || VectorQ[dim, (IntegerQ[#] && # > 0)&])

SDistribution/: Mean[SDistribution[g_, h_, alpha_, x0_, f0_], opts___?OptionQ] := Module[{dom,x, opti, res}, dom = Domain[SDistribution[g, h, alpha, x0, f0], FilterOptions[Domain,opts]]; If[dom===$Failed, Return[$Failed], With[{opti = FilterOptions[NIntegrate, opts]}, res = NIntegrate[ x PDF[SDistribution[g, h, alpha, x0, f0], x, opts], {x, dom[[1,1]], dom[[1,2]]}, opti]]; Return[res]]] /; ParameterQ[SDistribution[g, h, alpha, x0, f0]]



LerchPhiCNCT[z_, s_, v_, opts___?OptionQ] := Module[{reslist, stat, ret, it, p, p1},    
   p = LerchPhiRelativeAccuracy /. {opts} /. Options[LerchPhiCNCT]; 
   p1=1.0 10^(-p);
   reslist = Lerch[z,s,v,p1]; stat = reslist[[1]]; ret = reslist[[2]]; it = reslist[[3]];
   Which[
   stat==0, Return[ret], 
   stat==1, Message[LerchPhiCNCT::dvrg,stat]; Return[$Failed], 
   stat==2, Message[LerchPhiCNCT::intvdef,stat]; Return[$Failed], 
   stat==3, Message[LerchPhiCNCT::anegsr,stat]; Return[$Failed], 
   stat==4, Message[LerchPhiCNCT::intaj,stat]; Return[$Failed], 
   stat==5, Message[LerchPhiCNCT::underom,stat]; Return[$Failed], 
   stat==6, Message[LerchPhiCNCT::cnvg,stat]; Return[$Failed]
]] /; Head[z]=!=Complex && Head[s]=!=Complex && Head[v]=!=Complex

LerchPhiCNCT::dvrg = LerchDvrg1
LerchPhiCNCT::intvdef = LerchIntvdef2
LerchPhiCNCT::anegsr = LerchAnegsr3
LerchPhiCNCT::intaj = LerchIntaj4
LerchPhiCNCT::underom = LerchUnderom5
LerchPhiCNCT::cnvg = LerchCnvg6



(* Functions needed for calculations of quantiles (not exported). *)

LerchPhiClassical1[z_, s_, v_, sw1_, opts___?OptionQ] := Module[{reslist, stat, ret, it, p, p1}, 
   Which[sw1==1, p = LerchPhiRelativeAccuracy /. {opts} /. Options[Quantile], sw1==2, p = LerchPhiRelativeAccuracy /. {opts} /. Options[Random], sw1==3, p = LerchPhiRelativeAccuracy /. {opts} /. Options[RandomArray], sw1==4, p = LerchPhiRelativeAccuracy /. {opts} /. Options[Domain]];
   p1=1.0 10^(-p);
   reslist = Lerch[z,s,v,p1]; stat = reslist[[1]]; ret = reslist[[2]]; it = reslist[[3]];
   Which[
   sw1==1,
   Which[
   stat==0, Return[ret], 
   stat==1, Message[Quantile::dvrg,stat]; Return[$Failed], 
   stat==2, Message[Quantile::intvdef,stat]; Return[$Failed], 
   stat==3, Message[Quantile::anegsr,stat]; Return[$Failed], 
   stat==4, Message[Quantile::intaj,stat]; Return[$Failed], 
   stat==5, Message[Quantile::underom,stat]; Return[$Failed], 
   stat==6, Message[Quantile::cnvg,stat]; Return[$Failed]
   ],
   sw1==2,
   Which[
   stat==0, Return[ret], 
   stat==1, Message[Random::dvrg,stat]; Return[$Failed], 
   stat==2, Message[Random::intvdef,stat]; Return[$Failed], 
   stat==3, Message[Random::anegsr,stat]; Return[$Failed], 
   stat==4, Message[Random::intaj,stat]; Return[$Failed], 
   stat==5, Message[Random::underom,stat]; Return[$Failed], 
   stat==6, Message[Random::cnvg,stat]; Return[$Failed]
   ],
   sw1==3,
   Which[
   stat==0, Return[ret], 
   stat==1, Message[RandomArray::dvrg,stat]; Return[$Failed], 
   stat==2, Message[RandomArray::intvdef,stat]; Return[$Failed], 
   stat==3, Message[RandomArray::anegsr,stat]; Return[$Failed], 
   stat==4, Message[RandomArray::intaj,stat]; Return[$Failed], 
   stat==5, Message[RandomArray::underom,stat]; Return[$Failed], 
   stat==6, Message[RandomArray::cnvg,stat]; Return[$Failed]
   ],
   sw1==4,
   Which[
   stat==0, Return[ret], 
   stat==1, Message[Domain::dvrg,stat]; Return[$Failed], 
   stat==2, Message[Domain::intvdef,stat]; Return[$Failed], 
   stat==3, Message[Domain::anegsr,stat]; Return[$Failed], 
   stat==4, Message[Domain::intaj,stat]; Return[$Failed], 
   stat==5, Message[Domain::underom,stat]; Return[$Failed], 
   stat==6, Message[Domain::cnvg,stat]; Return[$Failed]
   ]]
] /; Head[z]=!=Complex && Head[s]=!=Complex && Head[v]=!=Complex

Quantile::dvrg = LerchDvrg1
Quantile::intvdef = LerchIntvdef2
Quantile::anegsr = LerchAnegsr3
Quantile::intaj = LerchIntaj4
Quantile::underom = LerchUnderom5
Quantile::cnvg = LerchCnvg6
Random::dvrg = LerchDvrg1
Random::intvdef = LerchIntvdef2
Random::anegsr = LerchAnegsr3
Random::intaj = LerchIntaj4
Random::underom = LerchUnderom5
Random::cnvg = LerchCnvg6
RandomArray::dvrg = LerchDvrg1
RandomArray::intvdef = LerchIntvdef2
RandomArray::anegsr = LerchAnegsr3
RandomArray::intaj = LerchIntaj4
RandomArray::underom = LerchUnderom5
RandomArray::cnvg = LerchCnvg6
Domain::dvrg = LerchDvrg1
Domain::intvdef = LerchIntvdef2
Domain::anegsr = LerchAnegsr3
Domain::intaj = LerchIntaj4
Domain::underom = LerchUnderom5
Domain::cnvg = LerchCnvg6



(* This version of the quantile calculation uses CNCT implementation of Lerch's Phi from loaded binary lerch.exe *)

SDistributionQuantile[{g_, h_, alpha_, x0_, f0_}, f_?NumberQ, sw_, opts___?OptionQ] := Module[{k, result, wasOn, lerch1, lerch2},    
wasOn = Head[Power::"infy"] =!= $Off; 
If[wasOn, Off[Power::"infy"]];  
k = (g 10-1.0 10)/(h 10-g 10);
result = 
Which[
(* nongeneric case g==1 *) 
IntegerPart[k]==k&&k==0,         
x0 + 1/alpha * (Log[f/f0] + 1/(h-1) * Log[ (1-f0^(h-1)) / (1-f^(h-1)) ] ), 
(* nongeneric case g > 1 *)  
IntegerPart[k]==k&&k>0,
x0 + 1/alpha * (Log[f/f0] + Sum[(f^(i*(h-g)-g+1) - f0^(i*(h-g)-g+1))/(i*(h-g)-g+1),{i,0,k-1}] + Sum[(f^(i*(h-g)-g+1) - f0^(i*(h-g)-g+1))/(i*(h-g)-g+1),{i,k+1,Infinity}]),
(* generic case *)
IntegerPart[k]!=k,
lerch1 = LerchPhiClassical1[f^(h-g),1.0,(1-g)/(h-g),sw,opts]; 
If[lerch1===$Failed, Return[$Failed]];
lerch2 = LerchPhiClassical1[f0^(h-g),1.0,(1-g)/(h-g),sw,opts]; 
If[lerch2===$Failed, Return[$Failed]];
x0 + 1/alpha/(h-g) * (f^(1-g)*lerch1 - f0^(1-g) * lerch2)];
If[wasOn, On[Power::"infy"]];
Return[result]] /; (sw==1 && (MethodPhi/.{opts}/.Options[Quantile])===LerchPhiCNCT) || (sw==2 && (MethodPhi/.{opts}/.Options[Random])===LerchPhiCNCT) || (sw==3 && (MethodPhi/.{opts}/.Options[RandomArray])===LerchPhiCNCT) || (sw==4 && (MethodPhi/.{opts}/.Options[Domain])===LerchPhiCNCT)



(* This version of the quantile calculation uses Mathematica implementation LerchPhi; for v<0 and integer s, I use an identity to transform to v>0  *)

SDistributionQuantile[{g_, h_, alpha_, x0_, f0_}, f_?NumberQ, sw_, opts___?OptionQ] := Module[{k, result, wasOn, lerch1, lerch2},    
wasOn = Head[Power::"infy"] =!= $Off; 
If[wasOn, Off[Power::"infy"]];  
k = (g 10-1.0 10)/(h 10-g 10);
result = 
Which[
(* nongeneric case g==1 *) 
IntegerPart[k]==k&&k==0,         
x0 + 1/alpha * (Log[f/f0] + 1/(h-1) * Log[ (1-f0^(h-1)) / (1-f^(h-1)) ] ), 
(* nongeneric case g > 1 *)  
IntegerPart[k]==k&&k>0,
x0 + 1/alpha * (Log[f/f0] + Sum[(f^(i*(h-g)-g+1) - f0^(i*(h-g)-g+1))/(i*(h-g)-g+1),{i,0,k-1}] + Sum[(f^(i*(h-g)-g+1) - f0^(i*(h-g)-g+1))/(i*(h-g)-g+1),{i,k+1,Infinity}]),
(* generic case *)
IntegerPart[k]!=k,
lerch1 = LerchPhiBuiltin[f^(h-g),1.0,(1-g)/(h-g),FilterOptions[LerchPhi,opts]]; 
If[lerch1===$Failed, Return[$Failed]];
lerch2 = LerchPhiBuiltin[f0^(h-g),1.0,(1-g)/(h-g),FilterOptions[LerchPhi,opts]]; 
If[lerch2===$Failed, Return[$Failed]];
x0 + 1/alpha/(h-g) * (f^(1-g)*lerch1 - f0^(1-g) * lerch2)];
If[wasOn, On[Power::"infy"]];
Return[result]] /; (sw==1 && (MethodPhi/.{opts}/.Options[Quantile])===LerchPhi) || (sw==2 && (MethodPhi/.{opts}/.Options[Random])===LerchPhi) || (sw==3 && (MethodPhi/.{opts}/.Options[RandomArray])===LerchPhi) || (sw==4 && (MethodPhi/.{opts}/.Options[Domain])===LerchPhi)

(* Use functional substitution like
/.LerchPhi[z_,s_,a_] :>NSum[z^k/(k+a)^s,{k,0,Infinity}]
*)

LerchPhiBuiltin[z_, s_, v_, opts___?OptionQ] := LerchPhi[z, s, v, opts] /; (v>0)
LerchPhiBuiltin[z_, s_, v_, opts___?OptionQ] := Module[{mm, i}, mm = -Floor[v]; z^mm LerchPhi[z, s, v + mm, opts] + Sum[z^i/(v + i)^s, {i, 0, mm - 1}]] /; (v<0)



(* Functions needed for calculation of cumulative distribution function (not exported). *)

SDistributionDF[par_List, x_List, sw_, opts___?OptionQ] := Module[{reslist, stat, ret, lprec, ztol, rtol, atol, meth, iter, jac, mm, steps, lprec1, atol1, rtol1}, 

   Which[
   sw==1,
   lprec = LerchPhiRelativeAccuracy /. {opts} /. Options[PDF];
   ztol = ZeroinIntervalUncertainty /. {opts} /. Options[PDF];
   atol = CVODEAbsoluteTolerance /. {opts} /. Options[PDF];
   rtol = CVODERelativeTolerance /. {opts} /. Options[PDF];
   meth = CVODEIntegrationMethod /. {opts} /. Options[PDF];
   iter = CVODEIterationType /. {opts} /. Options[PDF];
   jac = CVODEJacobian /. {opts} /. Options[PDF];
   steps = CVODEMaxSteps /. {opts} /. Options[PDF],
   sw==2,
   lprec = LerchPhiRelativeAccuracy /. {opts} /. Options[CDF];
   ztol = ZeroinIntervalUncertainty /. {opts} /. Options[CDF];
   atol = CVODEAbsoluteTolerance /. {opts} /. Options[CDF];
   rtol = CVODERelativeTolerance /. {opts} /. Options[CDF];
   meth = CVODEIntegrationMethod /. {opts} /. Options[CDF];
   iter = CVODEIterationType /. {opts} /. Options[CDF];
   jac = CVODEJacobian /. {opts} /. Options[CDF];
   steps = CVODEMaxSteps /. {opts} /. Options[CDF]];
   
   Which[meth === BDF, mm = 1000, meth === ADAMS, mm = 1100, True, mm = 1000];
   Which[iter === NEWTON, mm += 0, iter = FUNCTIONAL, mm += 10, True, mm += 0];
   Which[jac === ANALYTICAL, mm += 0, jac === NUMERICAL, mm += 1, True, mm += 0];

   lprec1=1.0 10^(-lprec); 
   atol1 = 1.0 10^(-atol);
   If[rtol==0,rtol1=0.0,rtol1=1.0 10^(-rtol)];
   reslist = SdCDF[par, x, lprec1, ztol, rtol1, atol1, mm, steps]; stat = reslist[[1]]; ret = reslist[[2]];
   
   If[stat==0, Return[ret]];
   
   Which[stat==1, Which[sw==1, Message[PDF::cvderr1, stat-9], sw==2, Message[CDF::cvderr1, stat-9]]; Return[$Failed], stat==2, Which[sw==1, Message[PDF::cvderr2, stat-9], sw==2, Message[CDF::cvderr2, stat-9]]; Return[$Failed], stat==3, Which[sw==1, Message[PDF::cvderr3, stat-9], sw==2, Message[CDF::cvderr3, stat-9]]; Return[$Failed], stat==4, Which[sw==1, Message[PDF::cvderr4, stat-9], sw==2, Message[CDF::cvderr4, stat-9]]; Return[$Failed], stat==5, Which[sw==1, Message[PDF::cvderr5, stat-9], sw==2, Message[CDF::cvderr5, stat-9]]; Return[$Failed], stat==6, Which[sw==1, Message[PDF::cvderr6, stat-9], sw==2, Message[CDF::cvderr6, stat-9]]; Return[$Failed], stat==7, Which[sw==1, Message[PDF::cvderr7, stat-9], sw==2, Message[CDF::cvderr7, stat-9]]; Return[$Failed], stat==8, Which[sw==1, Message[PDF::cvderr8, stat-9], sw==2, Message[CDF::cvderr8, stat-9]]; Return[$Failed]];
   
   Which[sw==1,
   Which[
   stat==9, Message[PDF::cvdneg]; Return[$Failed],
   stat==10, Message[PDF::dvrg,stat-9]; Return[$Failed], 
   stat==11, Message[PDF::intvdef,stat-9]; Return[$Failed], 
   stat==12, Message[PDF::anegsr,stat-9]; Return[$Failed], 
   stat==13, Message[PDF::intaj,stat-9]; Return[$Failed], 
   stat==14, Message[PDF::underom,stat-9]; Return[$Failed], 
   stat==15, Message[PDF::cnvg,stat-9]; Return[$Failed]
   ],
   sw==2,
   Which[
   stat==9, Message[CDF::cvdneg]; Return[$Failed],
   stat==10, Message[CDF::dvrg,stat-9]; Return[$Failed], 
   stat==11, Message[CDF::intvdef,stat-9]; Return[$Failed], 
   stat==12, Message[CDF::anegsr,stat-9]; Return[$Failed], 
   stat==13, Message[CDF::intaj,stat-9]; Return[$Failed], 
   stat==14, Message[CDF::underom,stat-9]; Return[$Failed], 
   stat==15, Message[CDF::cnvg,stat-9]; Return[$Failed]
   ]]

] /; ParameterQ[Apply[SDistribution,par]]

PDF::cvderr1 = SdCDFCvderr8
PDF::cvderr2 = SdCDFCvderr7
PDF::cvderr3 = SdCDFCvderr6
PDF::cvderr4 = SdCDFCvderr5
PDF::cvderr5 = SdCDFCvderr4
PDF::cvderr6 = SdCDFCvderr3
PDF::cvderr7 = SdCDFCvderr2
PDF::cvderr8 = SdCDFCvderr1
PDF::cvdneg = SdCDFCvdneg
PDF::dvrg = LerchDvrg1
PDF::intvdef = LerchIntvdef2
PDF::anegsr = LerchAnegsr3
PDF::intaj = LerchIntaj4
PDF::underom = LerchUnderom5
PDF::cnvg = LerchCnvg6
CDF::cvderr1 = SdCDFCvderr8
CDF::cvderr2 = SdCDFCvderr7
CDF::cvderr3 = SdCDFCvderr6
CDF::cvderr4 = SdCDFCvderr5
CDF::cvderr5 = SdCDFCvderr4
CDF::cvderr6 = SdCDFCvderr3
CDF::cvderr7 = SdCDFCvderr2
CDF::cvderr8 = SdCDFCvderr1
CDF::cvdneg = SdCDFCvdneg
CDF::dvrg = LerchDvrg1
CDF::intvdef = LerchIntvdef2
CDF::anegsr = LerchAnegsr3
CDF::intaj = LerchIntaj4
CDF::underom = LerchUnderom5
CDF::cnvg = LerchCnvg6



(* Functions needed for parameter estimation (exported). *)

SDistributionParameterEstimate[data_?VectorQ, opts___?OptionQ] :=
  
  Module[{emp, xx, ff, q75, f75, a, a1, x0, eq1, eq2, sol, f0, g, h, mmethod, 
      cums, freq, regrsol, x, n, k, ibin},
    
    AdaptiveBinCounts[data1_, in_] := 
      Module[{mn, mx, w, k, res, cat, pos, ibin, n},
        TakeMiddle[orig1_] := Take[orig1, {2, Length[orig1] - 1}];
        mx = Max[data1]; mn = Min[data1]; n = Length[data1]; 
        ibin = Round[n/in] + 1;
        w = {mn}; Do[w = {w, mn + k (mx - mn)/ibin}, {k, 1, ibin}];
        res = 
          NestWhile[(cat = TakeMiddle[RangeCounts[data1, #]]; 
                pos = Position[cat, 0]; 
                If[pos === {}, #, Drop[#, First[pos] + 1]]) &, Flatten[w], 
            UnsameQ, 2];
        Return[TakeMiddle[RangeCounts[data1, res]]]
        ];
    
    f0 = SampleQuantile /. {opts} /. Options[SDistributionParameterEstimate];
    If[f0 ² 0 || 1 ² f0, Message[SDistributionParameterEstimate::ffrac]; 
      Return[$Failed]];
    mmethod = 
      EstimateMethod /. {opts} /. Options[SDistributionParameterEstimate];
    If[mmethod =!= QuantileEquation && mmethod =!= AlgebraicRegression && 
        mmethod =!= ApproximateMLE, 
      Message[SDistributionParameterEstimate::wrme]; 
      mmethod = AlgebraicRegression];
    ibin = 
      NumberPointsInBins /. {opts} /. 
        Options[SDistributionParameterEstimate];
    If[! (IntegerQ[ibin] && Positive[ibin]), 
      Message[SDistributionParameterEstimate::wrib]; ibin = 6];
    
    a = 1/StandardDeviation[data]; x0 = Quantile[data, f0]; n = Length[data];
    
    Which[mmethod === QuantileEquation,
      
      emp = EmpiricalDF[data]; xx = emp[[1]]; ff = emp[[2]];
      q75 = Quantile[xx, 0.75]; f75 = ff[[Position[xx, q75][[1]]]][[1]];
      eq1 = Quantile[SDistribution[g, g + k, a, x0, f0], 0] == Min[data];
      eq2 = Quantile[SDistribution[g, g + k, a, x0, f0], f75] == q75;
      sol = 
        FindRoot[{eq1, eq2}, {g, 0.5, 0.6}, {k, 1.1, 2.1}, FilterOptions[FindRoot,opts],
          MaxIterations -> 100];
      Return[{g /. sol, g + k /. sol, a, x0, f0}],
      
      mmethod === AlgebraicRegression,
      
      freq = AdaptiveBinCounts[data, ibin]/n // N; 
      cums = CumulativeSums[freq];
      regrsol = 
        NonlinearRegress[Transpose[{cums, freq}], 
          a (x^g - x^(g + k)), {x}, {{g, 0.5, -100, 100}, {k, 0.01, 0, 100}}, FilterOptions[NonlinearRegress,opts],
          MaxIterations -> 1000];
      Return[{g, g + k, a, x0, f0} /. (BestFitParameters /. regrsol)],
      
      mmethod === ApproximateMLE,
      
      edf = Drop[EmpiricalDF[data][[2]], -1];
      eq1 = 1/(g + 1) + 1/(g + 1 + k) == -Apply[Plus, Log[edf]]/n;
      eq2 = (1 - n)/k/n + 
            1/(g + 1 + k) == -Apply[Plus, Log[edf]/(1 - edf^(-k))]/n;
      sol = 
        FindRoot[{eq1, eq2}, {g, 2, -10, 10}, {k, 0.01, 0, 10}, FilterOptions[FindRoot,opts],
          MaxIterations -> 100];
      Return[{g, h, a, x0, f0} /. h -> g + k /. sol],
      
      True, Return[$Failed]];
    
    ]

SDistributionParameterEstimate::ffrac =
"SampleQuantile is expected to be between 0 and 1."
SDistributionParameterEstimate::wrme = "EstimateMethod is not QuantileEquation, AlgebraicRegression, or ApproximateMLE; default value AlgebraicRegression is used instead."
SDistributionParameterEstimate::wrib = "NumberPointsInBins is not a positive integer; default value 6 is used instead."



SDistributionOptimize[data_?VectorQ, par:{{{_,_,_,_,_},_}...}, opts___?OptionQ] := Module[{symmode1, mode1, eps1, mit, reslist, stat, ret, statname, uid, lprec, ztol, atol, rtol, mm, steps, iter, meth, jac, passdt, delmin, delmax, delta, lowbnd, maxnf, maxit, anoise, rnoise, pivthr, addthr, xchthr, cnstol, npmin, layer, effort, varnt, scale, iprint, iout, lprec1, atol1, rtol1, par1, optimswitch, optsw, cpartmode, cfsqpndiffdelta, cfsqpobjeps, cfsqpobjrep, cfsqpglgeps},  

   delmin = DFOTrustLowerBound /. {opts} /. Options[SDistributionOptimize];
   delmax = DFOTrustUpperBound /. {opts} /. Options[SDistributionOptimize];
   delta = DFOTrustInitial /. {opts} /. Options[SDistributionOptimize];
   lowbnd = DFOObjectiveLowerBound /. {opts} /. Options[SDistributionOptimize];
   maxnf = DFOMaxFunctionEvaluations /. {opts} /. Options[SDistributionOptimize];
   maxit = DFOMaxIterations /. {opts} /. Options[SDistributionOptimize];
   anoise = DFOAbsoluteFunctionNoise /. {opts} /. Options[SDistributionOptimize];
   rnoise = DFORelativeFunctionNoise /. {opts} /. Options[SDistributionOptimize];
   pivthr = DFOPivotThreshold /. {opts} /. Options[SDistributionOptimize];
   addthr = DFOPivotThresholdGeometry /. {opts} /. Options[SDistributionOptimize];
   xchthr = DFOPivotFactor /. {opts} /. Options[SDistributionOptimize];
   npmin = DFOMinPointsInterpolation /. {opts} /. Options[SDistributionOptimize];
   layer = DFOMaxBaseDistanceFactor /. {opts} /. Options[SDistributionOptimize];
   effort = DFOInterpolationEffort /. {opts} /. Options[SDistributionOptimize];
   varnt = DFOQuadraticModel /. {opts} /. Options[SDistributionOptimize];
   scale = DFOScaling /. {opts} /. Options[SDistributionOptimize];
   iprint = -1; (*DFOPrintoutLevel /. {opts} /. Options[SDistributionOptimize];*)
   iout = 6; (*DFOOutputDevice /. {opts} /. Options[SDistributionOptimize];*)
   cfsqpobjeps = CFSQPObjectiveAbsoluteEpsilon /. {opts} /. Options[SDistributionOptimize];
   cfsqpobjrep = CFSQPObjectiveRelativeEpsilon /. {opts} /. Options[SDistributionOptimize];
   cfsqpglgeps = CFSQPKuhnTuckerFinalNorm /. {opts} /. Options[SDistributionOptimize];   
   eps1 = CFSQPNewtonDirectionFinalNorm /. {opts} /. Options[SDistributionOptimize];
   If[eps1 < 0, Message[SDistributionOptimize::negep]; Return[$Failed]];
   cfsqpndiffdelta = CFSQPFiniteDifferencePerturbation /. {opts} /. Options[SDistributionOptimize];
   If[cfsqpndiffdelta < 0, Message[SDistributionOptimize::negep]; Return[$Failed]];    
   symmode1 = CFSQPAlgorithm /. {opts} /. Options[SDistributionOptimize];
   cpartmode = CFSQPCheckConstraintsFirst /. {opts} /. Options[SDistributionOptimize];
	Which[cpartmode===True, mode1=200, cpartmode===False, mode1=100, True, Message[SDistributionOptimize::chker]; mode1=100];
   Which[symmode1===FSQPNL, mode1+=10, symmode1===FSQPAL, mode1+=0, True, Message[SDistributionOptimize::alg]; mode1+=10];
   mit = CFSQPMaxIterations /. {opts} /. Options[SDistributionOptimize];
   If[!NumericQ[mit] || !IntegerQ[mit] || !Positive[mit], Message[SDistributionOptimize::mitrng]; Return[$Failed]];
   optsw = OptimizationAlgorithm /. {opts} /. Options[SDistributionOptimize];
   statname = Statistic /. {opts} /. Options[SDistributionOptimize];
   optimswitch = Which[optsw===DFO, 0, optsw===CFSQP, 1, True, Message[SDIstributionOptimize::optmeth]; 0];
   uid = Which[
   statname===KolmogorovSmirnov, 0, 
   statname===Kuiper, 1, 
   statname===CramerVonMises, 2, 
   statname===AndersonDarling, 4, 
   statname===Watson, 3, 
   statname===RMS, 5, 
   True, Message[SDistributionOptimize::stat]; 4];   
   lprec = LerchPhiRelativeAccuracy /. {opts} /. Options[SDistributionOptimize];
   ztol = ZeroinIntervalUncertainty /. {opts} /. Options[SDistributionOptimize];
   atol = CVODEAbsoluteTolerance /. {opts} /. Options[SDistributionOptimize];
   rtol = CVODERelativeTolerance /. {opts} /. Options[SDistributionOptimize];
   meth = CVODEIntegrationMethod /. {opts} /. Options[SDistributionOptimize];
   iter = CVODEIterationType /. {opts} /. Options[SDistributionOptimize];
   jac = CVODEJacobian /. {opts} /. Options[SDistributionOptimize];
   steps = CVODEMaxSteps /. {opts} /. Options[SDistributionOptimize];   
   Which[meth === BDF, mm = 1000, meth === ADAMS, mm = 1100, True, mm = 1000];
   Which[iter === NEWTON, mm += 0, iter = FUNCTIONAL, m += 10, True, mm += 0];
   Which[jac === ANALYTICAL, mm += 0, jac === NUMERICAL, m += 1, True, mm += 0];

   lprec1 = 1.0 10^(-lprec);
   atol1 = 1.0 10^(-atol);
   If[rtol==0,rtol1=0.0,rtol1=1.0 10^(-rtol)];   
   par1 = Map[Flatten,par];   
   reslist = SdOptim[uid, Sort[data], par1, delmin, delmax, delta, lowbnd, maxnf, maxit, anoise, rnoise, pivthr, addthr, xchthr, npmin, layer, effort, varnt, scale, iprint, iout, lprec1, ztol, rtol1, atol1, mm, steps, eps1, cfsqpobjeps, cfsqpobjrep, cfsqpglgeps, cfsqpndiffdelta, mode1, mit, optimswitch]; 
   If[reslist===$Failed, Return[$Failed]];
   stat = reslist[[1]]; ret = {reslist[[2]], reslist[[3]]};
   
   If[stat == 25 || stat == 28, Return[ret]];   
   Which[
   stat==1, Message[SDistributionOptimize::cvderr1, stat-9]; Return[$Failed], 
   stat==2, Message[SDistributionOptimize::cvderr2, stat-9]; Return[$Failed], 
   stat==3, Message[SDistributionOptimize::cvderr3, stat-9]; Return[$Failed], 
   stat==4, Message[SDistributionOptimize::cvderr4, stat-9]; Return[$Failed], 
   stat==5, Message[SDistributionOptimize::cvderr5, stat-9]; Return[$Failed], 
   stat==6, Message[SDistributionOptimize::cvderr6, stat-9]; Return[$Failed], 
   stat==7, Message[SDistributionOptimize::cvderr7, stat-9]; Return[$Failed], 
   stat==8, Message[SDistributionOptimize::cvderr8, stat-9]; Return[$Failed],
   stat==9, Message[SDistributionOptimize::cvdneg]; Return[$Failed],
   stat==10, Message[SDistributionOptimize::dvrg, stat-9]; Return[$Failed], 
   stat==11, Message[SDistributionOptimize::intvdef, stat-9]; Return[$Failed], 
   stat==12, Message[SDistributionOptimize::anegsr, stat-9]; Return[$Failed], 
   stat==13, Message[SDistributionOptimize::intaj, stat-9]; Return[$Failed], 
   stat==14, Message[SDistributionOptimize::underom, stat-9]; Return[$Failed], 
   stat==15, Message[SDistributionOptimize::cnvg, stat-9]; Return[$Failed]
   ];
   Which[
   stat==16, Message[SDistributionOptimize::dfo1, stat-15-10]; Return[$Failed], 
   stat==17, Message[SDistributionOptimize::dfo2, stat-15-10]; Return[$Failed],
   stat==18, Message[SDistributionOptimize::dfo3, stat-15-10]; Return[$Failed],
   stat==19, Message[SDistributionOptimize::dfo4, stat-15-10]; Return[$Failed],
   stat==20, Message[SDistributionOptimize::dfo5, stat-15-10]; Return[$Failed],
   stat==21, Message[SDistributionOptimize::dfo6, stat-15-10]; Return[$Failed],
   stat==22, Message[SDistributionOptimize::dfo7, stat-15-10]; Return[$Failed],
   stat==23, Message[SDistributionOptimize::dfo8, stat-15-10]; Return[$Failed],
   stat==24, Message[SDistributionOptimize::dfo9, stat-15-10]; Return[$Failed],
   stat==26, Message[SDistributionOptimize::dfo11, stat-15-10]; Return[$Failed],
   stat==27, Message[SDistributionOptimize::dfo12, stat-15-10]; Return[$Failed]
   ];   
   Which[
   stat==29, Message[SDistributionOptimize::cfsqp1, stat-15-13]; Return[$Failed],
   stat==30, Message[SDistributionOptimize::cfsqp2, stat-15-13]; Return[$Failed],
   stat==31, Message[SDistributionOptimize::cfsqp3, stat-15-13]; Return[$Failed],
   stat==32, Message[SDistributionOptimize::cfsqp4, stat-15-13]; Return[$Failed],
   stat==33, Message[SDistributionOptimize::cfsqp5, stat-15-13]; Return[$Failed],
   stat==34, Message[SDistributionOptimize::cfsqp6, stat-15-13]; Return[$Failed],
   stat==35, Message[SDistributionOptimize::cfsqp7, stat-15-13]; Return[$Failed],
   stat==36, Message[SDistributionOptimize::cfsqp8, stat-15-13]; Return[$Failed],
   stat==37, Message[SDistributionOptimize::cfsqp9, stat-15-13]; Return[$Failed]
   ];
] (*/; DomainQ[Apply[SDistribution,par], Min[data]] *)

SDistributionOptimize::negep = 
"This value is expected to be positive."
SDistributionOptimize::stat =
"Warning: option Statistic is not specified correctly, using default value AndersonDarling."
SDistributionOptimize::alg = 
"Warning: option CFSQPAlgorithm is not specified correctly, using default value FSQPNL."
SDIstributionOptimize::optmeth =
"Warning: option OptimizationMethod is not specified correctly, using default value DFO."
SDistributionOptimize::mitrng =
"Value for option CFSQPMaxIterations is expected to be a positive integer."
SDistributionOptimize::chker =
"Warning: value for option CFSQPCheckConstraintsFirst is neither True or False, using default value False." 
SDistributionOptimize::cvderr1 = SdCDFCvderr8
SDistributionOptimize::cvderr2 = SdCDFCvderr7
SDistributionOptimize::cvderr3 = SdCDFCvderr6
SDistributionOptimize::cvderr4 = SdCDFCvderr5
SDistributionOptimize::cvderr5 = SdCDFCvderr4
SDistributionOptimize::cvderr6 = SdCDFCvderr3
SDistributionOptimize::cvderr7 = SdCDFCvderr2
SDistributionOptimize::cvderr8 = SdCDFCvderr1
SDistributionOptimize::cvdneg = SdCDFCvdneg
SDistributionOptimize::dvrg = LerchDvrg1
SDistributionOptimize::intvdef = LerchIntvdef2
SDistributionOptimize::anegsr = LerchAnegsr3
SDistributionOptimize::intaj = LerchIntaj4
SDistributionOptimize::underom = LerchUnderom5
SDistributionOptimize::cnvg = LerchCnvg6
SDistributionOptimize::dfo1 = SDOptimdfo9m
SDistributionOptimize::dfo2 = SDOptimdfo8m
SDistributionOptimize::dfo3 = SDOptimdfo7m
SDistributionOptimize::dfo4 = SDOptimdfo6m
SDistributionOptimize::dfo5 = SDOptimdfo5m
SDistributionOptimize::dfo6 = SDOptimdfo4m
SDistributionOptimize::dfo7 = SDOptimdfo3m
SDistributionOptimize::dfo8 = SDOptimdfo2m
SDistributionOptimize::dfo9 = SDOptimdfo1m
SDistributionOptimize::dfo11 = SDOptimdfo1p
SDistributionOptimize::dfo12 = SDOptimdfo2p
SDistributionOptimize::cfsqp1 = SDOptimcfsqp1
SDistributionOptimize::cfsqp2 = SDOptimcfsqp2
SDistributionOptimize::cfsqp3 = SDOptimcfsqp3
SDistributionOptimize::cfsqp4 = SDOptimcfsqp4
SDistributionOptimize::cfsqp5 = SDOptimcfsqp5
SDistributionOptimize::cfsqp6 = SDOptimcfsqp6
SDistributionOptimize::cfsqp7 = SDOptimcfsqp7
SDistributionOptimize::cfsqp8 = SDOptimcfsqp8
SDistributionOptimize::cfsqp9 = SDOptimcfsqp9



(* Miscelanea. *)

EmpiricalDF[data_?VectorQ, opts___?OptionQ] := Module[{edfv, currpoint, oldpoint, j, k, tieflag, edf, edf2, data1, data2, len}, 
   mode = EmpiricalDFDefinition /. {opts} /. Options[EmpiricalDF];
   cnst = EmpiricalDFDefinitionConstant /. {opts} /. Options[EmpiricalDF];
   If[ cnst<0 || cnst>1, Message[EmpiricalDF::cnst,cnst]; Return[$Failed]];
   If[ mode =!= KeepTies && mode =!= EliminateTies, Message[EmpiricalDF::tiesm,mode]; mode = EliminateTies];
   data1=Sort[data]; len=Dimensions[data1][[1]];
   Which[
      mode===EliminateTies,
    edf=Table[0.0,{len}]; edfv=1.0/len; edf[[1]]=edfv; tieflag=False; k=1; oldpoint=data1[[1]]; 
   For[j=2, j<=len, j++,   
   currpoint=data1[[j]]; edfv+=1.0/len;    
   If[((currpoint==oldpoint)&&(!tieflag)), (*k--;*) tieflag=True];   
   If[((currpoint!=oldpoint)&&(tieflag)), tieflag=False; (*k++;*) data1[[k]]=oldpoint; edf[[k]]=edfv-1.0/len]; 
   If[((currpoint==oldpoint)&&(j==len)), tieflag=False; (*k++;*) data1[[k]]=oldpoint; edf[[k]]=edfv]; 
   If[currpoint!=oldpoint, k++; data1[[k]]=currpoint; edf[[k]]=edfv];  
      oldpoint=currpoint;
     ]; (* For *)
   data2=Table[data1[[j]], {j,1,k}]; edf2=Table[edf[[j]], {j,1,k}]; 
   Return[{data2, edf2}],
      mode===KeepTies,
   Return[ {data1,Table[(j - cnst)/(len - 2 cnst + 1) // N, {j, 1, len}]}]
   ]
]

EmpiricalDF::cnst =
"Definition constant `1` is expected to be between 0 and 1."
EmpiricalDF::tiesm =
"Warning: Illegal option value `1`; using default value EliminateTies."



GoodnessOfFitStatistic[data_List, dist_, opts___?OptionQ] := Module[
   {statname, stat, listedf, dat, edf, cdf, n},   
   (* listedf = EmpiricalDF[data];
   dat = listedf[[1]]; edf = listedf[[2]]; *)
   n = Length[data];
   dat = Sort[data]; edf = Table[i, {i,1,n}] / n //N;
   cdf = CDF[dist,dat,FilterOptions[CDF,opts]];
   If[cdf===$Failed, Return[$Failed], Scan[If[#===$Failed, Return[$Failed]]&, cdf]];
   statname = Statistic /. {opts} /. Options[GoodnessOfFitStatistic];
   stat = Which[
   statname===KolmogorovSmirnov, Return[Max[Max[edf-cdf],Max[cdf-edf+1.0/n]]], 
   statname===Kuiper, Return[Max[edf-cdf] + Max[cdf-edf+1.0/n]], 
   statname===CramerVonMises, Return[Apply[Plus, (cdf-edf+1.0/2/n)^2]+1.0/12/n], 
   statname===AndersonDarling, Return[-n-Apply[Plus, (2 edf - 1.0/n) (Log[cdf] + Log[1.0-Reverse[cdf]])]], 
   statname===Watson, Return[Apply[Plus, (cdf-edf+1.0/2/n)^2] + 1.0/12/n - n (Apply[Plus, cdf]/n-1.0/2)^2], 
   statname===RMS, Return[Sqrt[Apply[Plus, (cdf-edf)^2]]], 
   True, Message[GoodnessOfFitStatistic::stat]; Return[-n-Apply[Plus, (2 edf - 1.0/n) (Log[cdf] + Log[1.0-Reverse[cdf]])]]];
]

GoodnessOfFitStatistic::stat =
"Warning: Statistic is not specified correctly, using default value AndersonDarling."



SampleQuantileCI[data1_?VectorQ, u_?QuantileQ, opts___?OptionQ] := Module[{qU1, qU2, n, u1, u2, u1s1, u1s2, u2s1, u2s2, eps1, eps2, cov, alpha, data}, 
   cov = ConfidenceLevel /. {opts} /. Options[SampleQuantileCI]; If[cov <= 0 || 1 <= cov, Message[SampleQuantileCI::clev]; Return[$Failed]]; alpha=1-cov;
   data = Sort[data1];
   n=Length[data];   
   u1s1=Quantile[BetaDistribution[(n+1)*u,(n+1)*(1-u)],alpha/2];  
   u1s2=Quantile[BetaDistribution[(n+1)*u,(n+1)*(1-u)],alpha/1.9];  
   u2s1=Quantile[BetaDistribution[(n+1)*u,(n+1)*(1-u)],1-alpha/2];  
   u2s2=Quantile[BetaDistribution[(n+1)*u,(n+1)*(1-u)],1-alpha/2.1];  
u1=u1/.FindRoot[BetaRegularized[u,(n+1)*u1,(n+1)*(1-u1)]==1-alpha/2,{u1,{u1s1,u1s2}}];  
u2=u2/.FindRoot[BetaRegularized[u,(n+1)*u2,(n+1)*(1-u2)]==alpha/2,{u2,{u2s1,u2s2}}]; eps1=(n+1)*u1-Floor[(n+1)*u1]; eps2=(n+1)*u2-Floor[(n+1)*u2];  
   qU1=(1-eps1)*data[[Floor[(n+1)*u1]]]+eps1*data[[Floor[(n+1)*u1]+1]];     
   qU2=(1-eps2)*data[[Floor[(n+1)*u2]]]+eps2*data[[Floor[(n+1)*u2]+1]]; 
   Return[{qU1,qU2}]
   ]

SampleQuantileCI::clev =
"ConfidenceLevel is expected to be between 0 and 1."



End[]



Protect[SDistribution];

(*, Quantile, Random, RandomArray, CDF, PDF, Domain, DomainQ, ParameterQ, Mean, Variance, StandardDeviation, Skewness, Kurtosis*)


EndPackage[]
