S distribution

This package provides Mathematica functions to calculate standard statistical functionals and parameter estimation (including  minimum distance) for the S distribution.

Author:  Sergej V. Aksenov

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

Last updated 3 May 2002

Installation and initialization

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

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

a second in the preferences directory:

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

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

Important note: 

this package requires three MathLink binaries lerch.exe, cdf.exe, and optim.exe to be present in the same directory as SDistribution.m.  
For some instructions on how to compile and build these binaries for your particular computer platform please refer to file README.

If the binaries cannot be found, or it is elected on purpose, the package can be used for calculation of the S distribution with standard Mathematica functions LerchPhi and NDSolve.  In this case the following options have to be turned on MethodODE->NDSolve and MethodPhi->LerchPhi.

The  SDistribution` package can  be loaded with

<<SDistribution`SDistribution`

Introduction

S distribution is defined by the following differential equation for its CDF

(1)         (dF/dx) = α ((F^g) - 
            (F^h))

where its parameters are subject to the following inequalities

α > 0,   g < h

ODE (1) can be used to solve for CDF (right hand side will give PDF) and when inverted for quantiles X.  Quantiles can also be calculated with Lerch's transcendent.

There have been four methods to estimate parameters of the S distribution:
* Algebraic regression
* Solving quantile equation
* Approximate maximum likelihood (ML) method
* alpha and x0 are inverse of standard deviation and sample quantile, respectively.

These are realized in function SDistributionParameterEstimate.

The main purpose of this package is to use Minimum distance estimation that is realized in function SDistributionOptimize.

Examples

Example 1

In Example 1 we calculate standard statistical functionals of the S distribution with provided MathLink binaries.

This defines the symbol SDistribution..

?SDistribution

We can examine if parameters of the S distribution are valid:

?ParameterQ

This is a well-defined S distribution..

ParameterQ[SDistribution[0.5, 1.6, 1.0, 0.0, 0.5]]

The following distributions are given ill-defined parameters..

ParameterQ[SDistribution[2.5, 1.6, 1.0, 0.0, 0.5]]

ParameterQ[SDistribution[0.5, 1.6, -1.0, 0.0, 0.5]]

ParameterQ[SDistribution[0.5, 1.6, 1.0, 0.0, 0.0]]

We can check what is the domain of the specified S distribution and if a given point is in this domain:

?Domain

This S distribution has a finite left tail..

Domain[SDistribution[0.5, 1.6, 1.0, 0.0, 0.5]]

This S distribution has an infinite left tail..

Domain[SDistribution[1.5, 1.6, 1.0, 0.0, 0.5]]

?DomainQ

Not all points are described by an S distribution with a finite left tail..

DomainQ[SDistribution[0.5, 1.6, 1.0, 0.0, 0.5], -1.0]

DomainQ[SDistribution[0.5, 1.6, 1.0, 0.0, 0.5], -10.0]

We can calculate extreme and intermediate quantiles:

?Quantile

Quantile[SDistribution[0.5, 1.6, 1.0, 0.0, 0.5], 0]

Quantile[SDistribution[0.5, 1.6, 1.0, 0.0, 0.5], 1]

Quantile[SDistribution[0.5, 1.6, 1.0, 0.0, 0.5], 0.1]

We can generate random numbers..

?Random

SeedRandom[879012]

Random[SDistribution[0.5, 1.6, 1.0, 0.0, 0.5]]

?RandomArray

SeedRandom[109283]

RandomArray[SDistribution[0.5, 1.6, 1.0, 0.0, 0.5], 10]

We can calculate PDF and CDF..

?PDF

Plot[
        PDF[SDistribution[0.5, 1.6, 1.0, 0.0, 0.5], x], {x, 
          Min[Domain[SDistribution[0.5, 1.6, 1.0, 0.0, 0.5]]], 
          10}];

?CDF

Plot[
        CDF[SDistribution[0.5, 1.6, 1.0, 0.0, 0.5], x], {x, 
          Min[Domain[SDistribution[0.5, 1.6, 1.0, 0.0, 0.5]]], 
          10}];

Example 2

In Example 2 we estimate parameters of the S distribution from univariate data.

Generate a sample from an S distribution.

SeedRandom[543967]

sddata = 
        RandomArray[SDistribution[0.4, 2.5, 1.0, 0.0, 0.5], 
          100];

Calculate estimates of S distribution parameters using algebraic regression method.

?SDistributionParameterEstimate

SDistributionParameterEstimate[sddata]

Since x0 is estimated as a sample quantile (median in this case), we can calculate the nonparametric 0.95 confidence interval..

Quantile[sddata, 0.5]

?SampleQuantileCI

SampleQuantileCI[sddata, 0.5]

Calculate goodness-of-fit statistic (Anderson-Darling) for the data and estimated S distribution..

?GoodnessOfFitStatistic

GoodnessOfFitStatistic[sddata, 
      Apply[SDistribution, {1.2850136115202087`, 1.596873554197161`, 
          1.4967298242205929`, -0.007797492451661947`, 0.5`}]]

Check if estimate is a valid set of parameters given the data..

DomainQ[
      Apply[SDistribution, {1.2850136115202087`, 1.596873554197161`, 
          1.4967298242205929`, -0.007797492451661947`, 0.5`}], 
      Min[sddata]]

Estimate S distribution parameters using minimum distance method (goodness-of-fit statistic is Anderson-Darling).

?SDistributionOptimize

SDistributionOptimize[
      sddata, {{{1.2850136115202087`, 1.596873554197161`, 
            1.4967298242205929`, -0.007797492451661947`, 0.5`}, 
          0.37963005358654256`}}]

Let's examine the quality of the fit visually.
For this we need to display an empirical CDF from the data:

?EmpiricalDF

DisplayTogether[ListPlot[Transpose[EmpiricalDF[sddata]]], 
        Plot[CDF[
            SDistribution[0.6504640686411545`, 2.0753916349212975`, 
              1.59822914366122`, 0.0027307996039558358`, 0.5`], x], {x, 
            Min[sddata], Max[sddata]}]];

Example 3

In Example 3 we calculate standard statistical functionals of the S distribution with the built-in Mathematica functions NDSolve and LerchPhi instead of the MathLink binaries.

Domain[SDistribution[0.5, 1.6, 1.0, 0.0, 0.5], 
      MethodPhi -> LerchPhi]

Quantile[SDistribution[0.5, 1.6, 1.0, 0.0, 0.5], 0]

Quantile[SDistribution[0.5, 1.6, 1.0, 0.0, 0.5], 0.1]

SeedRandom[879012]

Random[SDistribution[0.5, 1.6, 1.0, 0.0, 0.5]]

Plot[
        PDF[SDistribution[0.5, 1.6, 1.0, 0.0, 0.5], x, 
            MethodODE -> NDSolve] - 
          PDF[SDistribution[0.5, 1.6, 1.0, 0.0, 0.5], x], {x, 
          Min[Domain[SDistribution[0.5, 1.6, 1.0, 0.0, 0.5]]], 
          10}];

Plot[
        CDF[SDistribution[0.5, 1.6, 1.0, 0.0, 0.5], x, 
            MethodODE -> NDSolve] - 
          CDF[SDistribution[0.5, 1.6, 1.0, 0.0, 0.5], x], {x, 
          Min[Domain[SDistribution[0.5, 1.6, 1.0, 0.0, 0.5]]], 
          10}];

SeedRandom[543967]

sddata = 
        RandomArray[SDistribution[0.4, 2.5, 1.0, 0.0, 0.5], 
          100];

FindMinimum[
      GoodnessOfFitStatistic[sddata, 
        SDistribution[x, y, z, w, 0.5]], {x, {0.5, 0.6}}, {y, {2.5, 
          2.6}}, {z, {1.0, 1.1}}, {w, {0.1, 0.2}}]

SDistribution options

The following lists what options are available for each of the functions in the package.  
Options are then given below in alphabetical order and explained with examples; descriptions are sometimes adapted from the related external libraries' user's guides.

Options[Quantile]

Options[Random]

Options[RandomArray]

Options[PDF]

Options[CDF]

Options[DomainQ]

Options[Domain]

Options[SDistributionParameterEstimate]

Options[SampleQuantileCI]

Options[GoodnessOfFitStatistic]

Options[SDistributionOptimize]

Options[EmpiricalDF]

Options[LerchPhiCNCT]

Note: 

Examples of using some options below assume that data in Example 2 has been generated and is stored in sddata.

CFSQPAlgorithm

The following explains option CFSQPAlgorithm

?CFSQPAlgorithm

The following explains possible values for option CFSQPAlgorithm.

?FSQPNL

?FSQPAL

Default value is

CFSQPAlgorithm /. Options[SDistributionOptimize]

This option is included in mode in CFSQP code.  Please refer to CFSQP documentation for more information.

Example below shows that selecting FSQP-AL algorithm leads to problems with convergence.

SDistributionOptimize[
      sddata, {{{1.2850136115202087`, 1.596873554197161`, 
            1.4967298242205929`, -0.007797492451661947`, 0.5`}, 
          0.37963005358654256`}}, CFSQPAlgorithm -> FSQPAL]

CFSQPCheckConstraintsFirst

The following explains option CFSQPCheckConstraintsFirst

?CFSQPCheckConstraintsFirst

Possible values for option CFSQPCheckConstraintsFirst are True and False.

Default value is

CFSQPCheckConstraintsFirst /. Options[SDistributionOptimize]

This option is included in mode in CFSQP code.  Please refer to CFSQP documentation for more information.

Example below shows that selecting True does not seem to make any difference.

SDistributionOptimize[
      sddata, {{{1.2850136115202087`, 1.596873554197161`, 
            1.4967298242205929`, -0.007797492451661947`, 0.5`}, 
          0.37963005358654256`}}, 
      CFSQPCheckConstraintsFirst -> True]

CFSQPFiniteDifferencePerturbation

The following explains option CFSQPFiniteDifferencePerturbation

?CFSQPFiniteDifferencePerturbation

Possible values are nonnegative real numbers.

Default value is

CFSQPFiniteDifferencePerturbation /. 
      Options[SDistributionOptimize]

This option corresponds to udelta in CFSQP code.  The actual perturbation size is defined by Sign[x_i]*Max[udelta, Sqrt[macheps]*Max[1,Abs[x_i]] for each component of variable vector x_i, where macheps is the machine epsilon.  It is recommended to set it to zero if there is no specific reason to set it otherwise.  Please refer to CFSQP documentation for more information.

Example below shows that selecting a value other than 0.0 leads to problems with convergence.

SDistributionOptimize[
      sddata, {{{1.2850136115202087`, 1.596873554197161`, 
            1.4967298242205929`, -0.007797492451661947`, 0.5`}, 
          0.37963005358654256`}}, 
      CFSQPFiniteDifferencePerturbation -> 
        1.0 10^(-4)]

CFSQPKuhnTuckerFinalNorm

The following explains option CFSQPKuhnTuckerFinalNorm

?CFSQPKuhnTuckerFinalNorm

Possible values are positive real numbers.

Default value is

CFSQPKuhnTuckerFinalNorm /. Options[SDistributionOptimize]

This option corresponds to gLgeps in CFSQP code.  Set to a negative value (e.g., -1.0) not to use.  Please refer to CFSQP documentation for more information.

Example below shows that when the default stopping criterion (option CFSQPNewtonDirectionFinalNorm) is set to a unrealistically low value, convergence is governed by the Kuhn-Tucker norm.

SDistributionOptimize[
      sddata, {{{1.2850136115202087`, 1.596873554197161`, 
            1.4967298242205929`, -0.007797492451661947`, 0.5`}, 
          0.37963005358654256`}}, 
      CFSQPNewtonDirectionFinalNorm -> 1.0 10^(-14), 
      CFSQPKuhnTuckerFinalNorm -> 1.0 10^(-1)]

CFSQPMaxIterations

The following explains option CFSQPMaxIterations.

?CFSQPMaxIterations

Possible values are positive integers.

Default value is

CFSQPMaxIterations /. Options[SDistributionOptimize]

This option corresponds to miter in CFSQP code.    Please refer to CFSQP documentation for more information.

Example below shows that when the default stopping criterion (option CFSQPNewtonDirectionFinalNorm) is set to a low value, convergence is not achieved within the default number of iterations.  Allowing more iterations helps.

SDistributionOptimize[
      sddata, {{{1.2850136115202087`, 1.596873554197161`, 
            1.4967298242205929`, -0.007797492451661947`, 0.5`}, 
          0.37963005358654256`}}, 
      CFSQPNewtonDirectionFinalNorm -> 1.0 10^(-8)]

SDistributionOptimize[
      sddata, {{{1.2850136115202087`, 1.596873554197161`, 
            1.4967298242205929`, -0.007797492451661947`, 0.5`}, 
          0.37963005358654256`}}, 
      CFSQPNewtonDirectionFinalNorm -> 1.0 10^(-8), 
      CFSQPMaxIterations -> 1500]

CFSQPNewtonDirectionFinalNorm

The following explains option CFSQPNewtonDirectionFinalNorm.

?CFSQPNewtonDirectionFinalNorm

Possible values are positive real numbers.

Default value is

CFSQPNewtonDirectionFinalNorm /. 
      Options[SDistributionOptimize]

This option corresponds to eps in CFSQP code.    Please refer to CFSQP documentation for more information.

Example below shows that lowering its value, one obtains a more precise minimum of the objective function.

SDistributionOptimize[
      sddata, {{{1.2850136115202087`, 1.596873554197161`, 
            1.4967298242205929`, -0.007797492451661947`, 0.5`}, 
          0.37963005358654256`}}, 
      CFSQPNewtonDirectionFinalNorm -> 1.0 10^(-6)]

CFSQPObjectiveAbsoluteEpsilon

The following explains option CFSQPObjectiveAbsoluteEpsilon.

?CFSQPObjectiveAbsoluteEpsilon

Possible values are positive real numbers.

Default value is

CFSQPObjectiveAbsoluteEpsilon /. 
      Options[SDistributionOptimize]

This option corresponds to objeps in CFSQP code.  Set to a negative value (e.g., -1.0) not to use.  Please refer to CFSQP documentation for more information.

Example below shows that when the default stopping criterion (option CFSQPNewtonDirectionFinalNorm) is set to a unrealistically low value, convergence is governed by the objeps.

SDistributionOptimize[
      sddata, {{{1.2850136115202087`, 1.596873554197161`, 
            1.4967298242205929`, -0.007797492451661947`, 0.5`}, 
          0.37963005358654256`}}, 
      CFSQPNewtonDirectionFinalNorm -> 1.0 10^(-14), 
      CFSQPObjectiveAbsoluteEpsilon -> 1.0 10^(-3)]

CFSQPObjectiveRelativeEpsilon

The following explains option CFSQPObjectiveAbsoluteEpsilon.

?CFSQPObjectiveRelativeEpsilon

Possible values are positive real numbers.

Default value is

CFSQPObjectiveRelativeEpsilon /. 
      Options[SDistributionOptimize]

This option corresponds to objrep in CFSQP code.  Set to a negative value (e.g., -1.0) not to use.  Please refer to CFSQP documentation for more information.

Example below shows that when the default stopping criterion (option CFSQPNewtonDirectionFinalNorm) is set to a unrealistically low value, convergence is governed by the objrep.

SDistributionOptimize[
      sddata, {{{1.2850136115202087`, 1.596873554197161`, 
            1.4967298242205929`, -0.007797492451661947`, 0.5`}, 
          0.37963005358654256`}}, 
      CFSQPNewtonDirectionFinalNorm -> 1.0 10^(-14), 
      CFSQPObjectiveRelativeEpsilon -> 1.0 10^(-3)]

ConfidenceLevel

The following explains the standard Mathematica option ConfidenceLevel.

?ConfidenceLevel

Possible values are in the interval (0,1)

Default value is

ConfidenceLevel /. Options[SampleQuantileCI]

Example below shows that the nonparametric interval for sample median is wider when more confidence is demanded than 0.95.

SampleQuantileCI[sddata, 0.5]

SampleQuantileCI[sddata, 0.5, ConfidenceLevel -> 0.99]

CVODEAbsoluteTolerance

The following explains option CVODEAbsoluteTolerance.

?CVODEAbsoluteTolerance

Possible values are negative real numbers or zero, so that abstol = 
        (10^(-<option value>)).

Default value is

CVODEAbsoluteTolerance /. Options[PDF]

This option corresponds to abstol in CVODE code.  The estimated local error |(Subscript[e, i])| of the component (Subscript[y, i]) is less than reltol | (Subscript[y, i]) | +abstol.  reltol is set by the option CVODERelativeTolerance.  If abstol is set to 0, only reltol is used.  Please refer to CVODE documentation for more information.

The following example shows that solution of the S distribution ODE near the 0th quantile is sensitive to specification of absolute accuracy:

PDF[SDistribution[0.5, 1.6, 1.0, 0.0, 0.5], -1.7]

PDF[SDistribution[0.5, 1.6, 1.0, 0.0, 0.5], -1.7, 
      CVODEAbsoluteTolerance -> 6]

CVODEIntegrationMethod

The following explains option CVODEIntegrationMethod.

?CVODEIntegrationMethod

Possible values are:

?BDF

?ADAMS

Default value is

CVODEIntegrationMethod /. Options[PDF]

This option corresponds to lmm in CVODE code.  Select BDF for stiff ODE (which is apparently the case of S distribution ODE, at least in the tails of the distribution); ADAMS for nonstiff problems.  Please refer to CVODE documentation for more information.

The following example shows that solution of the S distribution ODE near the 0th quantile is sensitive to specification of solution method:

PDF[SDistribution[0.5, 1.6, 1.0, 0.0, 0.5], -1.707]

PDF[SDistribution[0.5, 1.6, 1.0, 0.0, 0.5], -1.707, 
      CVODEIntegrationMethod -> ADAMS]

CVODEIterationType

The following explains option CVODEIterationType.

?CVODEIterationType

Possible values are:

?NEWTON

?FUNCTIONAL

Default value is

CVODEIterationType /. Options[PDF]

This option corresponds to iter in CVODE code.  Select NEWTON for stiff ODE (which is apparently the case of S distribution ODE, at least in the tails of the distribution); FUNCTIONAL for nonstiff problems.  NEWTON uses the dense linear solver module of the CVODE.  Please refer to CVODE documentation for more information.

The following example shows that solution of the S distribution ODE is apparently insensitive to specification of iteration method:

PDF[SDistribution[0.5, 1.6, 1.0, 0.0, 0.5], 10.0]

PDF[SDistribution[0.5, 1.6, 1.0, 0.0, 0.5], 10.0, 
      CVODEIterationType -> FUNCTIONAL]

CVODEJacobian

The following explains option CVODEJacobian.

?CVODEJacobian

Possible values are:

?ANALYTICAL

?NUMERICAL

Default value is

CVODEJacobian /. Options[PDF]

This option corresponds to Jac in CVODE code.  Please refer to CVODE documentation for more information.

The following example shows that solution of the S distribution ODE is apparently insensitive to specification of jacobian:

PDF[SDistribution[0.5, 1.6, 1.0, 0.0, 0.5], -1.7]

PDF[SDistribution[0.5, 1.6, 1.0, 0.0, 0.5], -1.7, 
      CVODEJacobian -> NUMERICAL]

CVODEMaxSteps

The following explains option CVODEMaxSteps.

?CVODEMaxSteps

Possible values are positive integer numbers.

Default value is

CVODEMaxSteps /. Options[PDF]

This option corresponds to iopt[MXSTEP] in CVODE code.  Please refer to CVODE documentation for more information.

The following example shows that solution of the S distribution ODE near the 0th quantile is sensitive to specification of the maximum number of steps:

PDF[SDistribution[0.5, 1.6, 1.0, 0.0, 0.5], 25.0]

PDF[SDistribution[0.5, 1.6, 1.0, 0.0, 0.5], 25.0, 
      CVODEMaxSteps -> 100]

CVODERelativeTolerance

The following explains option CVODERelativeTolerance.

?CVODERelativeTolerance

Possible values are negative real numbers or zero, so that reltol = 
        (10^(-<option value>)).

Default value is

CVODERelativeTolerance /. Options[PDF]

This option corresponds to reltol in CVODE code.   The estimated local error |(Subscript[e, i])| of the component (Subscript[y, i]) is less than reltol | (Subscript[y, i]) | +abstol.  abstol is set by the option CVODEAbsoluteTolerance.  If reltol is set to 0, only abstol is used.  Please refer to CVODE documentation for more information.

The following example shows that solution of the S distribution ODE near the 0th quantile is sensitive to specification of relative accuracy:

PDF[SDistribution[0.5, 1.6, 1.0, 0.0, 0.5], -1.7]

PDF[SDistribution[0.5, 1.6, 1.0, 0.0, 0.5], -1.7, 
      CVODERelativeTolerance -> 6]

DFOAbsoluteFunctionNoise

The following explains option DFOAbsoluteFunctionNoise.

?DFOAbsoluteFunctionNoise

Possible values are positive real numbers.

Default value is

DFOAbsoluteFunctionNoise /. Options[SDistributionOptimize]

This option corresponds to anoise in DFO code.  Set to 0.0 not to use.  Please refer to DFO documentation for more information.

It is not clear what is the advantage of using this option with the S distribution.

SDistributionOptimize[
      sddata, {{{1.2850136115202087`, 1.596873554197161`, 
            1.4967298242205929`, -0.007797492451661947`, 0.5`}, 
          0.37963005358654256`}}, OptimizationAlgorithm -> DFO, 
      DFOAbsoluteFunctionNoise -> 1.0 10^(-3)]

DFOInterpolationEffort

The following explains option DFOInterpolationEffort.

?DFOInterpolationEffort

Possible values are integers 1-4.

Default value is

DFOInterpolationEffort /. Options[SDistributionOptimize]

This option corresponds to effort in the DFO code. This option specifies the level of effort the algorithm spends on the geometry.
Please refer to DFO documentation for more information.

It is not clear what is the advantage of using this option with the S distribution.  Sometimes it can save the situation.

SDistributionOptimize[
      sddata, {{{1.2850136115202087`, 1.596873554197161`, 
            1.4967298242205929`, -0.007797492451661947`, 0.5`}, 
          0.37963005358654256`}}, OptimizationAlgorithm -> DFO, 
      DFOInterpolationEffort -> 4]

DFOMaxBaseDistanceFactor

?DFOMaxBaseDistanceFactor

Possible values are positive real numbers.

Default value is

DFOMaxBaseDistanceFactor /. Options[SDistributionOptimize]

This option corresponds to layer in the DFO code.

Points that are DFOMaxBaseDistanceFactor times DFOTrustInitial away from the base are considered last.

Large DFOMaxBaseDistanceFactor allows models to be more complete but less up-to-date.  Low DFOMaxBaseDistanceFactor leads up-to-date but incomplete models (so information may be lost).  Values 2-10 are recommended.

Please refer to DFO documentation for more information.

Find example?

DFOMaxFunctionEvaluations

?DFOMaxFunctionEvaluations

Possible values are positive integers.

Default value is

DFOMaxFunctionEvaluations /. Options[SDistributionOptimize]

This option corresponds to maxnf in DFO code.  Note that single complete general quadratic interpolation requires (n+1)*(n+2)/2 function values, i.e. 15.

DFOMaxIterations

?DFOMaxIterations

Possible values are positive integers.

Default value is

DFOMaxIterations /. Options[SDistributionOptimize]

This value correpsonds to maxit in DFO code.  The real number of iterations is in fact limited by the number of function evaluaitons..

DFOMinPointsInterpolation

?DFOMinPointsInterpolation

Possible values are positive integers.

Default value is

DFOMinPointsInterpolation /. Options[SDistributionOptimize]

This option corresponds to npmin in DFO code.  The interpolation will include at least this number of points.  For example, if it is n+1 (5 for four S distribution parameters), then every interpolation will at least be fully linear, if at all possible.

DFOObjectiveLowerBound

?DFOObjectiveLowerBound

Possible values are real numbers.

Default value is

DFOObjectiveLowerBound /. Options[SDistributionOptimize]

This option corresponds to lowbnd in DFO code.  Useful if there is target value of objective function.

DFOPivotFactor

The following explains option DFOPivotFactor.

?DFOPivotFactor

Possible values are positive real numbers.

Default value is

DFOPivotFactor /. Options[SDistributionOptimize]

This option corresponds to xchthr in the DFO code.

This option specifies the factor of improvement in the pivot, when replacing an existing interpolation point by another one, with hopes to improve the geometry.  When geometry needs improvement, a point with with the low pivot is selected.  If the pivot produced by the new point is greater or equal to DFOPivotFactor times DFOTrustInitial, then the new point is accepted.

Large DFOPivotFactor restricts the number of point exchanges.  Low DFOPivotFactor leads to acceptance of points that do not improve geometry.

Please refer to DFO documentation for more information.

Find example?

SDistributionOptimize[sddata,{{{1.2850136115202087`, 1.596873554197161`, 
              1.4967298242205929`, -0.007797492451661947`, 0.5`}, 
            0.37963005358654256`}},OptimizationAlgorithm -> DFO,DFOPivotFactor->0.001]

DFOPivotThreshold

?DFOPivotThreshold

Possible values are positive real numbers.

Default value is

DFOPivotThreshold /. Options[SDistributionOptimize]

This option corresponds to pivthr in DFO code.  The lower the threshold, the larger the risk of ill-conditioned inaccurate models and the larger is the number of points that will be included in the interpolation.

DFOPivotThresholdGeometry

?DFOPivotThresholdGeometry

Possible values are positive real numbers.

Default value is

DFOPivotThresholdGeometry /. Options[SDistributionOptimize]

This option corresponds to addthr in DFO code.  This option works when a new point in the interpolation set does not produce any or good reduction; then it is still included if it produces a good pivot (determined by addthr).
Note that DFOPivotThresholdGeometry>=DFOPivotThreshold

DFOQuadraticModel

?DFOQuadraticModel

Possible values are 1 and 2.

Default value is

DFOQuadraticModel /. Options[SDistributionOptimize]

This option correspods to varnt in DFO code.  Effect of this option is not fully tested yet.  Generally Frobenius norm (1) is thought to work better.

DFORelativeFunctionNoise

The following explains option DFORelativeFunctionNoise.

?DFORelativeFunctionNoise

Possible values are positive real numbers.

Default value is

DFORelativeFunctionNoise /. Options[SDistributionOptimize]

This option corresponds to rnoise in DFO code.  Set to 0.0 not to use.  Please refer to DFO documentation for more information.

It is not clear what is the advantage of using this option with the S distribution.

SDistributionOptimize[
      sddata, {{{1.2850136115202087`, 1.596873554197161`, 
            1.4967298242205929`, -0.007797492451661947`, 0.5`}, 
          0.37963005358654256`}}, OptimizationAlgorithm -> DFO, 
      DFORelativeFunctionNoise -> 1.0 10^(-8)]

DFOScaling

?DFOScaling

Possible values are 0,1,2

Default value is

DFOScaling /. Options[SDistributionOptimize]

This option corresponds to scale in DFO code.  Scaling is recommended only if ncecssary.  Note that DFOTrustLowerBound is scaled accordingly, but DFOTrustInitial and DFOTrustUpperBound are not.

DFOTrustInitial

?DFOTrustInitial

Possible values are positive real numbers.

Default value is

DFOTrustInitial /. Options[SDistributionOptimize]

This value corresponds to delta in DFO code.

DFOTrustLowerBound

?DFOTrustLowerBound

Possible values are positive real numbers.

Default value is

DFOTrustLowerBound /. Options[SDistributionOptimize]

This option corresponds to delmin in DFO code.  This is the stopping criterion--whenever the trust region radius falls below this value the algorithm stops.  Values below 10^(-6) are not recommended.

DFOTrustUpperBound

?DFOTrustUpperBound

Possible values are positive real numbers.

Default value is

DFOTrustUpperBound /. Options[SDistributionOptimize]

This option corresponds to delmax in DFO code.  This bound helps prevent the trust region from growing too large and making too large a step, which makes old information irrelevant and new points have to be calculated.

EmpiricalDFDefinition

?EmpiricalDFDefinition

?EliminateTies

?KeepTies

EmpiricalDFDefinition /. Options[EmpiricalDF]

EmpiricalDFDefinitionConstant

?EmpiricalDFDefinitionConstant

EmpiricalDFDefinitionConstant /. Options[EmpiricalDF]

EstimateMethod

?EstimateMethod

Possible values are the following

?AlgebraicRegression

?QuantileEquation

?ApproximateMLE

Default value is

EstimateMethod /. Options[SDistributionParameterEstimate]

It seems that AlgebraicRegression is most useful value, except possibly left-skewed data.

SDistributionParameterEstimate[sddata]

SDistributionParameterEstimate[sddata, 
      EstimateMethod -> QuantileEquation]

SDistributionParameterEstimate[sddata, 
      EstimateMethod -> ApproximateMLE]

LerchPhiRelativeAccuracy

?LerchPhiRelativeAccuracy

Possible values are positive real numbers less than one.

Default value is

LerchPhiRelativeAccuracy /. Options[Quantile]

Value of this option is passed to external routine LERCHPHI.  It is probably not reasonable to obtain more than 14 digits of accuracy with CNCT.  Also, see implementation notes for some issues with the loss of precision for some parameter values.

Please note that this option is used in other functions (see list above), so it is important to set its value consistently throughout all computations.

MethodPhi

?MethodPhi

The values relevant to S distribution calculations are

?LerchPhi

?LerchPhiCNCT

This example shows that for near-degenerate S distributions (g and h almost on the degeneracy line) both Mathematica and CNCT implementations in standard precision give incorrect result.

The reason is the loss of precision for v<0 very close to integer.  Thus the remedy is to perform calculations of Lerch's Phi in extended precision.

(2 1.6 + 1)/(2 + 1)

Degeneracy constant k is near integer (S distribution is near-degenerate):

ScientificForm[(1.40000000000001 - 1)/(1.6 - 1.40000000000001), 
      16]

Quantile with the CNCT implementation of LerchPhi

Quantile[SDistribution[1.40000000000001, 1.6, 1.0, 0.0, 0.5], 
      0.01]

Quantile with Mathematica implementation of LerchPhi  (agrees with CNCT)

Quantile[SDistribution[1.40000000000001, 1.6, 1.0, 0.0, 0.5], 0.01, 
      MethodPhi -> LerchPhi]

Both are wrong however;  here's result using Mathematica implementation of Lerch Phi in quadruple precision:

Quantile[SDistribution[1.40000000000001`64, 1.6`64, 1.0, 0.0, 0.5], 
      0.01, MethodPhi -> LerchPhi]

For uses of this option with other SDistribution functions see Example 3..

MethodODE

?MethodODE

The values relevant to S distribution calculations are

?NDSolve

?CVODE

In the example below, we can see that results from CVODE and NDSolve are very different around the problematic zeroth quantile.

CDF[SDistribution[0.5, 1.6, 1.0, 0.0, 0.5], -1.7068]

CDF[SDistribution[0.5, 1.6, 1.0, 0.0, 0.5], -1.7068, 
      MethodODE -> NDSolve]

For uses of this option with other SDistribution functions see Example 3..

NumberPointsInBins

?NumberPointsInBins

Possible values are positive integers.

Default value is

NumberPointsInBins /. 
      Options[SDistributionParameterEstimate]

This option specifies the desired number of points per bin from which the number of bins to be used is calculated; the algorithm joins bins if they are empty...  The end result is a sequence of nonempty bins with possibly nonequal length.

This option is related to function SDistributionParameterEstimate with the option EstimateMethod->AlgebraicRegession (which see).

OptimizationAlgorithm

?OptimizationAlgorithm

Possible values are

?CFSQP

?DFO

Default value is

OptimizationAlgorithm /. Options[SDistributionOptimize]

Note that for CFSQP, only one initial point is necessary; but it has to be feasible (see Example 2).

For DFO, one point can be given and it doesn't have to be feasible (because the algorithm will reflect it with respect to bounds if the point is not feasible).

SDistributionOptimize[
      sddata, {{{1.2850136115202087`, 1.596873554197161`, 
            1.4967298242205929`, -0.007797492451661947`, 0.5`}, 
          0.37963005358654256`}}, 
      OptimizationAlgorithm -> DFO]

However, the result will be random (because the algorithm generates the second required point randomly).  (Evaluation below uses exactly the same starting point...)

SDistributionOptimize[
      sddata, {{{1.2850136115202087`, 1.596873554197161`, 
            1.4967298242205929`, -0.007797492451661947`, 0.5`}, 
          0.37963005358654256`}}, 
      OptimizationAlgorithm -> DFO]

Alternatively, two or more points can be given but they have to be feasible; the result is no longer random..

GoodnessOfFitStatistic[sddata, 
      SDistribution[1.4, 1.7, 1.1, 0.002, 0.5`]]

SDistributionOptimize[
      sddata, {{{1.2850136115202087`, 1.596873554197161`, 
            1.4967298242205929`, -0.007797492451661947`, 0.5`}, 
          0.37963005358654256`}, {{1.4, 1.7, 1.1, 0.002, 0.5`}, 
          0.37963005358654256`}}, 
      OptimizationAlgorithm -> DFO]

SampleQuantile

?SampleQuantile

Possible values are real numbers between 0 and 1.

Default value is

SampleQuantile /. Options[SDistributionParameterEstimate]

Statistic

?Statistic

Possible values are the folowing

?AndersonDarling

?KolmogorovSmirnov

?Kuiper

?CramerVonMises

?Watson

?RMS

Default value is

Statistic /. Options[GoodnessOfFitStatistic]

Note that this option is also used in SDistributionOptimize so it is important to use the same value consistently throughout all computations.

ZeroinIntervalUncertainty

?ZeroinIntervalUncertainty

Possible value are positive real numbers less than one.

Default value is

ZeroinIntervalUncertainty /. Options[SDistributionOptimize]

This option corresponds to tol in the Zeroin code.  It specifies uncertainty in the result x according to the formula 4*macheps*abs(x)+tol.

Please note that this option is also used in functions CDF and PDF, soo it is important to set its value consistently throughout all computations.

Implementation notes

Quantile etc. and Lerch's transcendent

One important implementation detail concerns calculation of S distribution quantiles.  These involve calculation of Lerch's transcendent (for general case of the quantile solution g>1).

General issue:

Calculation of k=(g-1)/(h-g) is subject of loss of precision if k is integer;  digits are lost when calculating differences in numerator and denominator.  Currently the hack is to calculate (g 10 - 10)/(h 10 - g 10).  Accurate k is important because it decides whether to use generic or nongeneric solutions.

Issues with Mathematica implementation:

1. For v<0 and odd s, LerchPhi calculates differently from the "classical" definition of Phi (because Mathematica uses "alternative" definition).  See discussion at http://library.wolfram.com/mathgroup/archive/1997/Oct/msg00118.html

LerchPhiCNCT[0.1, 1.0, -1.5]

LerchPhi[0.1, 1.0, -1.5]

Thus, in the code Phi is calculated with LerchPhi for v>0 and with transformation for v<0.

Sum[z^i/(v + i) /. z -> 0.1 /. v -> -1.5, {i, 0, 
          1}] + 0.1^2 LerchPhi[0.1, 1.0, 2 - 1.5]

2.  For v<0 very close to integer (i.e., for nearly degenerate S distributions) there is loss of precision in LerchPhi and the result is thus incorrect.

Sum[z^i/(v + i) /. z -> 0.1 /. 
          v -> -3.00000000000001, {i, 0, 3}] + 
      0.1^4 LerchPhi[0.1, 1.0, 4 - 3.00000000000001]

The remedy should be to request evaluation of Quantile and others in extended precision.

Sum[z^i/(v + i) /. z -> 0.1`64 /. 
          v -> -3.00000000000001`64, {i, 0, 3}] + 
      0.1`64^4 LerchPhi[0.1`64, 1.0`64, 
          4.0`64 - 3.00000000000001`64]

Issues with CNCT implementation:

The only problem is the above one for v<0 very close to integer, also loss of precision.  Using binary lerch.exe is not recommended in this case.  In future versions there will be lerch.exe with extended precision..

Note that calculation of S distribution quantiles with the CNCT implementation of Phi (the same as in lerch.exe) in the optim.exe is not subject to the above warnings (i.e., there is no loss of precision), since in optim.exe quantiles are calculated only for g<1 (i.e., v>0).

Troubleshooting

SDistributionOptimize errors

If CFSQP gives error 4 or 5 or DFO gives error -7 (min. over trust region failed, which is CFSQP related), try perturbing the starting point (happens sometimes with AndersonDarling); however, be careful to still have the point feasible.