Jeff Miller
Department of Psychology
University of Otago
Dunedin, New Zealand

DistFit: A Program to Fit Probability Distributions to Data Sets

DistFit Version 1.1 - Cupid Version 2.1- Documentation Version Jul 29, 2005

Contents

1  Overview

2  Detailed Introduction

2.1  Parameter Estimation Methods

Maximum Likelihood Estimation

For a given probability distribution, the likelihood (L) of a set of observed data values is

L = N Õ i=1  f(Xi)

where f is the distribution's probability density function (PDF), and the values X_i, i=1¼N, are the observed data values. For maximum likelihood estimation, the parameters of the distribution are chosen to be whatever parameter values maximize the value of L for the given data set. Thus, this estimation method uses all of the observed data values.

Minimum Chi-square Estimation

As is explained in section 5.3.4, the "chi-square goodness of fit test" is a standard way of assessing the fit of a predicted theoretical distribution to an observed data set. For this test, the observed data set is divided into K mutually exclusive and exhaustive bins, and the number of observed scores in each bin is simply counted. This tabulation of the data is called a "histogram." Then, this observed histogram is compared against a theoretically predicted histogram (i.e., bin counts) in order to compute a measure of the deviation between the predicted and observed distributions. This measure is called the observed chi-square value, c²_o, with larger values indicating larger discrepancies between the predicted and observed distributions. For minimum chi-square estimation, then, the parameter values of the distribution are chosen to be whichever ones produce the smallest value of c²_o. Thus, this estimation method uses a histogram of the observed data values (i.e., a count of the number of data values in each of a set of distinct bins).

Moment-based Estimation

With this estimation method, the distribution's parameter values are chosen to match the observed data (as closely as possible) with respect to the mean (first moment), variance (second moment), and possibly 3rd, 4th, etc central moments as well. Thus, this estimation method uses only these moment-based summaries of the observed data set.

Percentile Matching Estimation

With this estimation method, the distribution's parameter values are chosen to match the observed data (as closely as possible) with respect to the values at a certain specified set of percentiles of the distribution. Thus, this estimation method uses only a set of percentiles of the observed data set.

2.2  A Simple Example

C> DistFit@MLE100b

* A sample file illustrating the format of DistFit input files. 
* This file illustrates maximum-likelihood fitting with the data 
* contained in a separate file (MLE100b.Dat). 

* The comment following each parameter indicates what information 
* it conveys to the program. 

NDists 2               * Two distributions are to be fit to the data. 
Normal(0,1)          * The first distribution is the normal, 
Laplace(0,1)         * and the second distribution is the Laplace. 

* In the preceding two lines, 0 and 1 are guesses for the parameters. 
* For each distribution, the program adjusts these to the best-fitting 
* values, but it is important to choose reasonable guesses. 

TEXTOUTPUT   * This line instructs the program to write the *.txt output file. 
TABOUTPUT    * This line instructs the program to write the *.tab output file. 

MLEFIT MLE100b.Dat  * The fitting is to be done by MLE (Max likelihood). 
* The data to be fit are contained in the file MLE100b.dat 
* This is a plain ASCII file with one data value per line.

• The first line shows the best-fitting parameter values.
• The next line shows the obtained value of -2 Ln(Likelihood) corresponding to the best fit. By definition, the maximum likelihood value is the one that minimizes -2 Ln(Likelihood).
• The next set of lines, down to and including the ChiSquare line, report the results of a chi-square goodness of fit test comparing the observations against the fitted distribution. The column "Range" refers to ranges of observed values. The column "Observed" indicates how many observed values were in each range, and the column "Predicted" indicates how many values are predicted to be in each range by the distribution with the best-fitting parameter values. The ChiSquare value at the bottom gets larger as the fit is poorer; if the p value at the end of the line is less than 0.05, then it should probably be concluded that the observed values do not come from this theoretical distribution.
• The last line simply records the program's default setting to look for the best parameter estimates with a plain simplex search; this setting and other possible settings are explained in section 7.4.

Best-fitting parameters for distribution 1: Normal(0.05253,0.8675) 
-2 Ln(Likelihood):   127.6766 
Range    Observed      Predicted 
-4.0708 to  -1.0593:         16          9.997 
-1.0593 to  -0.6774:          4         10.008 
-0.6774 to  -0.4020:          7         10.009 
-0.4020 to  -0.1669:         13         10.002 
-0.1669 to   0.0525:          8          9.984 
0.0525 to   0.2719:         12          9.984 
0.2719 to   0.5071:         10         10.002 
0.5071 to   0.7825:         13         10.009 
0.7825 to   1.1644:          8         10.008 
1.1644 to   4.1758:          9          9.997 
ChiSquare =    11.2119 with 7 degrees of freedom (p=0.130) 
Search mode for distribution 1: Simplex 

Best-fitting parameters for distribution 2: LaPlace(0.1239,0.6844) 
-2 Ln(Likelihood):   131.3999 
Range    Observed      Predicted 
-12.0096 to  -0.9777:         17         10.000 
-0.9777 to  -0.5033:         10         10.000 
-0.5033 to  -0.2258:         10         10.000 
-0.2258 to  -0.0288:          8         10.000 
-0.0288 to   0.1239:          5         10.000 
0.1239 to   0.2766:         10         10.000 
0.2766 to   0.4735:         10         10.000 
0.4735 to   0.7510:         12         10.000 
0.7510 to   1.2255:          9         10.000 
1.2255 to  12.2574:          9         10.000 
ChiSquare =     8.4000 with 7 degrees of freedom (p=0.299) 
Search mode for distribution 2: Simplex

3  Installation & Test

1. Unzip the distribution file. It will be named something like DistFit.ZIP, except that the last few letters of DistFit will be changed to numbers to reflect the current version (e.g., DistFit11.ZIP is version 1.1). Note that subdirectories called Examples\DistFit\-In and Examples\DistFit\-Out.OK are created; these are for test purposes (see point after next).
2. Move the program file DistFit.EXE to a directory in your path. This program runs in a CMD window (also known as a DOS window), at least under Windows 2000 and XP.
3. If you like, you can run some quick tests of the program to make sure that everything is working properly on your machine and your version of the operating system. To do the quick tests, open a command window and change its current directory into the subdirectory Examples\DistFit\-In. Then run the batch file Go.Bat. The batch file should run a series of test runs, and it should produce output files as described later. When the batch file is done, check that every newly created file in the Examples\DistFit\-In subdirectory is identical to the file with the same name in the Examples\DistFit\-Out.OK subdirectory. If the files match, it is likely that everything is OK. If it fails, contact the author.

4  Starting the Program

C> DistFit@Exampl1 
C> DistFit@Exampl1a @Exampl1b 
C> DistFit@Exampl2 TabOutput
C> DistFit@Exampl1 MLEFit_MyDataFileName.Dat
C> DistFit@GotoExample1(My1stGotoLabel)

• Parameters may be specified on the command line itself, as in the third and fourth example DistFit commands shown above. In the third example command, for instance, DistFit reads parameters from Exampl2.Rsp and then sets one further parameter called TabOutput (described later).
• You can read program control parameters from two (or more) Rsp files, simply listing all of them on the command line. In the example on the second line, for instance, program control parameters are read from both Exampl1a.Rsp and @Exampl1b.Rsp.
• Sometimes control parameters in the Rsp file consist of two or more terms separated by white space, as you will see descriptions of the different control parameters. To specify such parameters on the command line, use an underscore instead of white space, as shown in the fourth command (MLEFit_MyDataFileName.Dat).
• You may want to use a single Rsp file to control several slightly different runs, with somewhat different parameter settings. When the Rsp file name is followed on the command line by an extra parameter in parentheses [e.g., @GotoExample1(My1stGotoLabel)], parameters are read from the indicated Rsp file starting at the point in the file indicated by the label in parentheses, rather than starting at the beginning of the Rsp file. The method of labelling a point in the Rsp file is explained further in section 5.5.

5  Structure of the Input File

• No parameters are case-sensitive.
• With some crucial exceptions, one parameter is specified on each line of the input file.
• The parameters can appear in (almost) any order.
• Most parameters are optional and have reasonable defaults.
• Parameters can be followed on the same line by comments. By default, the asterisk is used to start a comment, but this can be changed.
• Indenting (leading whitespace) is almost always optional (i.e., it is optional unless clearly noted otherwise). It is used in the examples to improve readability and indicate grouping.

5.1  Specifying the Distribution(s) to Fit

NDists  4            * This indicates that 4 distributions should be fit. 
Normal(0,1)        * Distribution 1 is normal, starting with mean=0 and sigma=1. 
Uniform(-10,10)    * Distribution 2 is uniform, starting from -10 to 10. 
t(15)              * Distribution 3 is a t distribution with 15 degrees of freedom. 
Normal(0,1)     FR * Distribution 4 is normal, starting with mean=0 and sigma=1.

5.2  Specifying the Estimation Method and the Data

5.2.1  Maximum Likelihood Estimation

MLEFit  MLE100b.dat   * Perform max likelihood fit to data in MLE100b.dat

MLEFit  4   * Perform max likelihood fit to following 4 data values 
0.123 
0.456 
0.777 
0.989

5.2.2  Minimum Chi-square Estimation

ChiSqFit  4   * Perform min chi-sq fit to data set with 4 bins as follows: 
0           * Bottom value of first bin 
10 4        * 10 is top value of 1st bin; the bin  0-10 has  4 data values. 
20 12       * 20 is top value of 2nd bin; the bin 10.0001-20 has 12 data values. 
30 12       * 30 is top value of 3rd bin; the bin 20.0001-30 has 12 data values. 
40 4        * 40 is top value of 4th bin; the bin 30.0001-40 has  4 data values.

5.2.3  Moment-based Estimation

MomentFIT 2  * Number of moments in this data set 
0.500     * The desired mean 
1.500     * The desired variance

5.2.4  Percentile Matching Estimation

PctileFit 5  * There are 5 percentile points to be fit. 
0.1 1       * The data value at the 10th percentile is 1 
0.3 2       * The data value at the 30th percentile is 2 
0.5 2.5     * The data value at the 50th percentile is 2.5 
0.7 3       * The data value at the 70th percentile is 3 
0.9 4       * The data value at the 90th percentile is 4

5.3  Controlling the Output

5.3.1  Choosing Output File(s)

5.3.2  Adjusting Number Widths and Numbers of Decimal Places

ParmWid  10 
ParmDP    5

5.3.3  Setting the Output File Name

DistFit@Test1(Goto1)

Outfile foo2

AppendOutfile MLEa

Outfile foo2 
AppendOutfile MLEa

5.3.4  Controlling Bins for the Chi-Square Goodness of Fit Test

ChiSquareGOFNBins 15       * Compute the chi-square(s) using 15 bins

ChisquareGOFBinsEqualInX   * Construct chi-square(s) bins equally wide in scores

ChisquareGOFBinsEqualInP   * Construct chi-square(s) bins to have equal probability

5.3.5  Adding a Title to the Output File

Title This is an example title for my output file.

5.4  Comments in the Input File

Comment !          * Starting with the NEXT LINE, comments are set off by ! 
Outfile MyOutName  ! Set the output file name 

5.5  Multiple Specifications in a Single Rsp File

TypeA            * Label for first set of specifications. 
...              * Here is the set of specifications for one run. 
End              * End of first set of specifications. 

TypeB            * Label for second set of specifications. 
...              * Here is the set of specifications for another run. 
Goto TypesBandC  * Go to the indicated label to process more commands. 
End              * End of second set of specifications. 

TypeC            * Label for third set of specifications. 
...              * Here is the set of specifications for another run. 
Goto TypesBandC  * Go to the indicated label to process more commands. 
End              * End of third set of specifications. 

TypesBandC       * Label that is the target for the above GOTO commands. 
...              * A set of commands to be used for both TypeB and TypeC. 
End

DistFitMultSpec(TypeA) 
DistFitMultSpec(TypeB)

5.6  Multiple Data Sets in a Single Run

5.6.1  Columnar Data for Maximum Likelihood Estimation

ColumnFormat            * Data is in a separate file in columnar format. 
MLEFit ColExMLE.Col     * ColExMLE.Col has 1 line per dataset.

220 180 800 540 670 320 280 600 940 710 
221 182 803 544 1010 833 675 326 287 608 949 710

5.6.2  Columnar Data for Minimum Chi-Square Estimation

ColumnFormat                   * Data is in a separate file in columnar format. 
BinEdges 200 250 300 350 400 450 500 550 600   * List of the high bin edge values 
MLEFit ColExChiSq.Col          * ColExChiSq.Col has 1 line per dataset.

1 10 13 16 22 13 11  9  5 
2 10 13 16 22 13 11  9  4

5.6.3  Columnar Data for Moment-Based Estimation

ColumnFormat            * Data is in a separate file in columnar format. 
MomentFit ColExMom.Col  * ColExMom.Col has 1 line per dataset.

240 2000 128000 
250 2100 134000

5.6.4  Columnar Data for Percentile-Based Estimation

ColumnFormat                      * Data is in a separate file in columnar format. 
Percentiles 0.1 0.3 0.5 0.7 0.9   * Data values are for these percentiles. 
PercentileFit ColExPct.Col        * ColExPct.Col has 1 line per dataset.

220  250  280  320 380 
220  250  285  330 410

6  Program Output

Sum of squared errors:     0.0016

7  Miscellaneous Advanced Options

7.1  Interrupting DistFit

7.2  Running DistFit at Low Priority

C> start /low DistFit@MyParms

7.3  ParmCodes

• The ParmCodes string should have exactly one character for each parameter in the distribution.
• Each character in the ParmCodes string should be either F, R, or I, with no other characters allowed.
• The order of the characters in the string should match the left-to-right order of the parameters in the distribution name.
• there is more complete information on ParmCodes setting in the documentation for CUPID.

Normal(0,1) FR

7.4  Search Modes

Simplex

This search mode uses the simplex algorithm of Rosenbrock (1960).

Grid search

This search mode tries all combinations of parameter values on a user-defined grid, possibly starting a new simplex search from each of the grid points.

Random search

This search mode tries randomly generated combinations of parameter values within user-defined bounds, possibly starting a new simplex search from each randomly generated combinations.

7.4.1  Specifying the Grid Search Mode

GridSearch 1* Use a grid search for search type 1. 
3 0.5 1.1     * Special line: See notes below at @@ 
100 500 10    * Parameter 1 should range from 100 to 500 in 10 steps. 
10 200 10    * Parameter 2 should range from 10 to 200 in 10 steps. 
1 400 10    * Parameter 3 should range from 1 to 400 in 10 steps. 

* @@ Notes for the above special line: 
* 3 indicates that there are three free parameters in this distribution. 
* 0.5 indicates that a simplex search should sometimes be started from a grid point. 
*     The alternatives are: 
*          0.0  Never start simplex search from a grid point. 
*          1.0  Always start simplex search from a grid point. 
* 1.1 indicates that the simplex search should start from a grid point if the error 
*     at that grid point is no more than CurrentBestError*1.1, where CurrentBestError 
*     is the best error found so far.

• The first number specifies the number of parameters to be varied in the search (more generally, the number of free parameters in the distribution being fitted to the data).
• The second number is used to specify whether the simplex search should be called at each of the grid points. There are three possible cases:
  1. If this value is less than or equal to 0, simplex is never called for any grid point.
  2. If this value is greater than or equal to 1, simplex is always called for every grid point.
  3. If this value is between 0 and 1, simplex is sometimes called, depending on the third number on this line.
• The third number, which will be called "ErrorFactor" is only used when the second number is between 0 and 1. In those cases, DistFit first evaluates the error associated with the current combination of parameters; call this error value "CurrentError". It then computes the ratio of CurrentError divided by ErrorFactor. If this ratio is less than the best error obtained previously, then DistFit  starts a simplex search at the current point; otherwise, it skips the simplex for this point. I usually use ErrorFactor values of around 1.1-1.5, with higher values doing more simplex searches (thus, taking more time, but also giving a better chance of finding the best overall result). Intuitively, what is going on here? Well, CurrentError measures the error at the current parameter combination. If that is fairly small, then it seems like we might be close to the right parameter combination so it would be worthwhile to start a simplex search from this point to refine it further. How small is "fairly small"? That is defined by comparing the current error to the best error score so far.

7.4.2  Specifying the Random Search Mode

RandomSearch 1
3 0.5 1.1 1000    * See notes below at @@ 
100 500           * Parameter 1 should range from 100 to 500. 
10 200           * Parameter 2 should range from 10 to 200. 
1 400           * Parameter 3 should range from 1 to 400. 

* @@ Notes for the above line: 
* 3 indicates that there are three free parameters in this distribution. 
* 0.5 indicates that a simplex search should sometimes be started from a grid point. 
*     The alternatives are: 
*          0.0  Never start simplex search from a grid point. 
*          1.0  Always start simplex search from a grid point. 
* 1.1 indicates that the simplex search should start from a grid point if the error 
*     at that grid point is no more than CurrentBestError*1.1, where CurrentBestError 
*     is the best error found so far. 
* 1000 indicates that 1000 starting points are to be randomly generated.

7.4.3  Special Output with Grid and Random Searches

Results: 620 searches 
 68 new bests found.

7.4.4  Simplex Search Parameters

MinStepSize

This option is used to control the minimum step size used by the simplex parameter search algorithm. For example:

MinStepSize 0.00001

The default is 1.0e-8.

SimplexStartingStep

This option is used to control the starting step size used by the simplex parameter search algorithm. For example:

SimplexStartingStep 0.01

The default is 0.2.

SimplexMaxIteration

This option is used to restrict the maximum number of points examined by the simplex parameter search algorithm. For example:

SimplexMaxIteration 10000

would allow up to 10000 points in the parameter space to be examined in a given simplex search. The default is 5000.

7.5  Controlling the Accuracy of Numerical Integration and Inversion

InverseprecisionX   0.001  * Compute inverses to within an X value of .001 (lenient). 
InverseprecisionP   0.001  * Compute inverses to within a P value of .001 (lenient).

7.6  Checking DistFit

8  Available Probability Distributions

8.1  Continuous Distributions

Beta(A,B)

The Beta distribution is defined over the interval from zero to one, and its shape is determined by its two parameters A and B. Its PDF is

f(x) = 1 b(A,B) xA-1 (1-x)B-1,    0 < x < 1

The mean is A/(A+B), and the variance is AB(A+B)^-2(A+B+1)^-1.

Cauchy(L,S)

This distribution is defined in terms of location and scale parameters L and S > 0, respectively. Its PDF is

f(x) = 1 pS é ë 1+ ì í î x-L S ü ý þ 2   ù û

ChiSquare(df)

This is a generalization of the distribution of the sum of df independent squared standard normals. Its parameter is df - a positive real number.

Chi(df)

This is the distribution of the positive square root of a ChiSquare random variable. Its parameter is df - a positive real number.

Cosine

For 0 £ x £ P/2, f(x)=cos(x) and F(x)=sin(x).

ExGaussian(m,s,l)

This is the distribution of the sum of independent Normal and Exponential random variables. Its parameters are the m and s of the Normal, and the rate l of the Exponential.

ExGaussMn(m,s,m_e)

This is a reparameterization of the ExGaussian. Its parameters are the m and s of the Normal, and the mean of the exponential, m_e.

ExGaussRat(m,s,r)

This is just a reparameterization of the ExGaussian. Its parameters are the m and s of the Normal, and the ratio, r, of the mean of the exponential to the sigma of the normal.

Exponential(l)

This distribution is well-known. By default, the parameter is the rate l; the mean is 1/l.

ExpSum(r1,r2)

This is the sum (convolution) of two exponentials with different rates. The two parameters are the two rates, which must be different enough to avoid numerical errors. For the convolution of exponentials with the same rates, of course, you should use the Gamma.

ExpSumT(r1,r2,Cutoff)

This is the sum (convolution) of two exponentials with different rates truncated at a given cutoff value. The first two parameters are the two rates, which must be different enough to avoid numerical errors; the third parameter is the upper truncation point. For the convolution of exponentials with the same rates, of course, you should use the Gamma.

ExpoNo(m,s)

I just invented this as an ad-hoc solution for a problem I was working on one time. I don't know whether it has ever been considered before or will ever be useful again, and I certainly don't know whether I gave it a reasonable name. In any case, it is a transformation of a normal random variable X. Specifically, it is the distribution of

Y = eX 1+eX

where X has a normal distribution with mean m and standard deviation s. The two parameters of this distribution are the m and s of the underlying normal X.

ExtremeVal(a,b)

Extreme-value Type I distribution (a.k.a. Fisher-Tippett distribution, Gumbel distribution, sometimes also called the double exponential distribution, to be confused with the Laplace distribution), with parameters a and b > 0. The CDF is

F(x) = exp[-e-(x-a)/b]

ExWald(m,s,a,l)

This is the distribution of the sum of two independent random variables: one from a three-parameter Wald distribution with parameters (m,s,a); and one from an exponential distribution with rate l. Schwarz (2001, 2002) describes the ex-Wald distribution in detail.

ExWaldMn(m,s,a,m_e)

This is a reparameterization of the ExWald. The first three parameters are the same as in the plain ExWald, and the fourth parameter is the mean of the exponential, m_e, instead of the rate.

ExWaldMSM(m,sd,m_e)

This is a further reparameterization of the ExWald. The first two parameters are the mean and standard deviation (not s!) of the Wald component, and the third parameter is the mean of the exponential, m_e. The Wald parameter a is set to 1.0.

F(dfNumer,dfDenom)

This is Fisher's distribution of the ratio of two independent normed Chi-square distributions, as commonly used in linear models (e.g., analysis of variance). The two integer parameters are the degrees of freedom of the numerator and denominator, respectively.

Gamma(N,l)

This is the distribution of the sum of N exponentials, each with rate l. In this distribution, N must be a positive integer. In the RNGamma distribution (see below), N is any positive real.

Geary(SampleSize)

The Geary statistic arises in testing to see whether a set of observations come from a normal distribution (D'Agostino, 1970).

GenErr(Mu,Scale,Shape)

This is the general (a.k.a. "generalized") error distribution (e.g., Evans, Hasting, & Peacock, 1993, p. 57), also known as Subbotin's distribution (e.g., Johnson, Kotz, & Balakrishnan, 1995, 2nd Ed., Vol 2, p. 195). The correct PDF is

f(x) = exp[-|x-Mu|Shape/(2·Scale)] Scale1/Shape ·2(1+1/Shape) ·G(1+1/Shape)

although both EHP and JKB give it with incorrect exponents of the Scale parameter in the denominator. This version uses the shape parameter denoted l by Evans et al. Note that the Laplace, normal, and uniform distributions are special cases of this distribution with this shape parameter equal to 1, 2, and approaching ¥, respectively. In practice, lots of combinations of parameter values give overflow errors, especially if the shape parameter is more than approximately 3.

HypTan(Scale)

This is the Hyperbolic Tangent distribution, whose PDF and CDF are

f(x) = 4·b [ebx + e-bx]2 F(x) = ebx - e-bx ebx + e-bx

where b is the scale parameter. This distribution arises as a model of psychometric functions (e.g., Strasburger, 2001).

Inverse Gaussian

See the Wald3 distribution.

Laplace(L,S)

Also known as the double exponential. In terms of location and scale parameters, L and S > 0, respectively, the PDF is

f(x) = 1 2 S e-|x-L|/S

Logistic(m,b)

This distribution is defined in terms of a location parameter m and a scale parameter b. The cumulative form of the distribution is

F(x) = 1 1 + e[(-(x-m))/(b)]

LogNormal(m,s)

This is the distribution of X such that ln(X) is normally distributed. The parameters are the m and s of the normal.

Naka-Rushton(Scale)

This is the distribution of X ³ 0 such that

f(x) = 2 ·x ·a2 (1+(a·x)2)2 F(x) = (a·x)2 1+(a·x)2

where a is the scale parameter. In the actual distribution, moments above the first do not exist; they do exist in DistFit's truncated version of the distribution, however.

NoncentralF(dfNumer,dfDenom,Noncentrality)

This is the distribution of the ratio of independent noncentral and central chi-squares, with the former in the numerator. It is most often used in the computation of power of the F test. The noncentrality parameter is defined in terms of the dfNumer normal random variables whose sum of squares is yields the chi-square in the numerator. Specifically,

l = dfNumer å i=1  Li2

where L_i is the expected value of the ith random variable contributing to this sum of squares.

Normal(m,s)

I'll bet you know this one already. Parameters are m and s, not s².

Pareto1(K,A)

This is a Pareto distribution of the first kind, as defined by Johnson, Kotz, and Balakrishnan (1994, vol 1, p 574), with PDF and CDF

f(x) = A ·KA ·x-(A+1) F(x) = 1 - æ è K A ö ø A

where K > 0, A > 0, and x ³ K. This is a model for income, where K is some minimum income and F(x) is the probability that a randomly selected income is less than or equal to x.

Quantal(Threshold)

This distribution is related to the Poisson. This is the distribution of X ³ 0 such that

F(x) = 1 - T-1 å t=0  xt t! e-x

This distribution arises as a model of psychometric functions in visual psychophysics (e.g., Gescheider, 1997, p. 85). The threshold parameter, T ³ 1, represents an observer's fixed threshold for the number of quanta of light that must be detected before saying "Yes, I saw the stimulus." Quanta are assumed to be emitted from the stimulus according to a Poisson distribution with parameter x. Then, F(x) is the psychometric function for the probability of saying "Yes" as a function of the mean number of quanta, x, emitted by the stimulus. Note that it makes no real sense to think of x as a random variable in this example, but the probability distribution provides a useful model anyway.

Quick(Scale,Shape)

This is the distribution of X ³ 0 with PDF and CDF

f(x) = 2-([x/(a)])b · æ è x a ö ø b   ·b·ln(2) x F(x) = 1 - 2-([x/(a)])b

where a is the scale parameter and b is the shape parameter. This distribution arises as a model of psychometric functions (e.g., Quick, 1974; Strasburger, 2001).

Rayleigh(s)

If Y₁ and Y₂ are independent normal random variables with mean 0 and standard deviation s, then X = Ö{Y₁²+Y₂²} has a Rayleigh distribution with scale parameter s. The PDF is

f(x) = e-x2/(2s2) x s2

RNGamma(RN,l)

See "Gamma". In this version, the shape parameter RN is a real number rather than an integer.

rPearson(SampleSize)

This is the sampling distribution of Pearson's r (correlation coefficient) under the null hypothesis that the true correlation is zero (and assuming the usual bivariate normality). The parameter is SampleSize, the number of pairs of observations across which the correlation is computed.

StudRng(df,NSamples)

Distribution of Studentized range statistic with df degrees of freedom for error and NSamples samples. Because both parameters are integers, automatic program-based estimation of these parameters is rarely successful.

t(df)

Student's t-distribution, with degrees of freedom equal to df.

Triangular(B,T)

In this distribution the density function has the shape of an equilateral triangle across some range. The parameters are the bottom (B) and the top of the range (T). The PDF is then:

f(x) = ì ï ï í ï ï î (x - B) ×Hp if B £ x £ B+T 2 (T - x) ×Hp if B+T 2 £ x £ T

where H_p is the height of the PDF at its peak, adjusted to so that the total area of the triangle is 1.0.

TriangularG(B,P,T)

In this (more general triangular) distribution, the density function has the shape of a not-necessarily-equilateral triangle across some range. The parameters are the bottom of the range (B), the point at which the triangle reaches its maximum (P), and the top of the range (T). The PDF is then:

f(x) = ì ï ï í ï ï î (x - B) ×Hp P-B if B £ x £ P (T - x) ×Hp T-P if P £ x £ T

where H_p is the height of the PDF at its peak, adjusted to so that the total area of the triangle is 1.0.

Uniform(B,T)

This is the distribution in which all values are equally likely within some range. The parameters are the bottom and the top of the range, B and T.

UniGap(T)

This is an equal-probability mixture of two uniform distributions, one extending from -T to 0 and the other extending from T to 2·T. It is "model 4" of Sternberg and Knoll (1973). The median is somewhat arbitrarily defined as T/2.

Wald3(m, s, a)

This is the general, three-parameter version of the Wald distribution. Specifically, assume a one-dimensional Wiener diffusion process starting at position 0 at time 0 and drifting with average rate m and variance s², and consider X to be the first passage time through position a. The PDF of X is

f(x) = a s Ö 2 px3 · exp é ë - (a - mx)2 2 s2 x ù û

where all three parameters and x must be positive. Note: By default, the sigma parameter is fixed in all parameter searches.

Weibull(Scale,Power,Origin)

As defined by Johnson & Kotz (1970, p. 250): "X has a Weibull distribution if there are values of the parameters c( > 0), a( > 0), and n₀ such that

Y = é ë (X-n0) a ù û c

has the exponential distribution with rate = 1". Here, the parameters c, a, and n₀ are referred to as the "scale," "power," and "origin" parameters, respectively.

8.2  Discrete Distributions

Binomial(N,p)

The distribution of the number of successes in N Bernoulli trials, with probability p of success on each trial.

Constant(C)

This is a degenerate distribution that always takes on the same value. Its parameter is that value. Perhaps surprisingly, it can be convenient to have this distribution available. Warning: For technical reasons, the constant distribution does not work well in many of the derived distributions discussed in the next section. Thus, it should be avoided whenever possible. For example, you should always use: LinearTrans(Gamma(2,.01),1,100) rather than the equivalent Convolution(Gamma(2,.01),Constant(100))

Geometric(P)

The distribution of the trial number of the first success in a sequence of Bernoulli trials, where P is the probability of success on each try.

List(filename)

This random variable allows you to define any discrete set of X values, each with its own arbitrary probability. The command List(FileName) tells DistFit to read the distribution from the indicated file. The first line in the file contains the number of X values in the distribution. After that, there should be one line for each X value, with the first number on the line being the X value itself and the second number being the probability of that X value. These values need not be sorted, and in fact the same X value can appear on several different lines, in which case the associated probabilities will be summed. (If X's do appear on more than one line, then the first line in the file should actually contain the number of X-containing lines rather than the number of distinct X's.)

Poisson(U)

X has a Poisson distribution with parameter U if

Pr (X=x) = e-UUx x! ,    x=0,1,2,¼, U > 0

The mean and variance both equal U.

UniformInt(Low,High)

This is the distribution of equally likely integer values between the two integer parameters, Low and High, inclusive.

8.3  Transformation Distributions

8.4  Derived Distributions

Convolution(RV1,RV2)

This is the distribution of a sum of independent random variables, RV1 and RV2, where RV1 and RV2 are each legal distributions in their own right. For example, Convolution(Normal(0,1),Uniform(0,1)) specifies the convolution of these normal and uniform distributions, and Convolution(Normal(0,1),Uniform(0,1),Gamma(3,0.01)) specifies the convolution of the three indicated distributions.

In general, to define a convolution, the user types something of the form: Convolution(BasisDist1(Parms),...,BasisDistK(Parms)) where Parms stands for the parameters associated with each of the distributions. There are K random variables summed together, and the distributions of these summed variables are simply listed, separated by commas.

DistFit is not very smart about convolutions. At this point, it only knows how to compute means, variances, and random numbers in an intelligent way. Everything else is computed using (recursive) numerical integration, which tends to be pretty slow. Also, DistFit does not "realize" that some convolutions result in a new distribution about which it already knows (e.g., convolution of two normals is normal). Thus, computations involving these convolutions proceed via numerical integration even though direct computation would be possible.

The current version can handle convolutions where all distributions are discrete, all are continuous, or some are is discrete and some continuous, but it cannot handle convolutions in which one or more distributions are mixed (i.e., partly discrete and partly continuous).

I would be very happy for suggestions on how to augment DistFit's handling of convolutions, especially those accompanied by Pascal code.

ConvolutionIID(N,RV)

This is just an easier way to specify a convolution when all N summed random variables have the same distribution, RV. ConvolutionIID(3,Uniform(0,1)) is the same as Convolution(Uniform(0,1),Uniform(0,1),Uniform(0,1))

Difference(RV1,RV2)

This is the distribution of the difference of two independent random variables, RV1 minus RV2, where RV1 and RV2 are each legal distributions in their own right. For example, Difference(Uniform(0,1),Uniform(0,1)) specifies a difference between two standard uniform distributions, which ranges from -1 to 1 (not uniformly). DistFit handles difference distributions dumbly, like convolutions. The current version can handle differences where both distributions are discrete, both are continuous, or one is discrete and one continuous, but it cannot handle differences in which one or both distributions are mixed (i.e., partly discrete and partly continuous).

Mixture(p1,RV1,p2,RV2,...,pk,RVk)

Mixtures are distributions formed by randomly selecting one of a number of random variables. For example, Mixture(0.5,Normal(0,1),0.5,Uniform(0,1)) defines a random variable that comes from a standard normal half the time and a standard uniform the other half of the time. In general, the format of this distribution is: Mixture(p₁,BasisDist1(Parms),p₂,BasisDist1(Parms),...,p_k,BasisDistK(Parms)) and the p_i's must sum to one (it is also legal to omit p_k).

InfMix(RV1,MixParm,RV2(Parms))

The InfMix distribution is an infinite mixture, formed when a parameter of one distribution is itself randomly distributed according to another distribution. For example, InfMix(Normal(0,5),1,Uniform(10,20)) defines a random variable that comes from a normal distribution with standard deviation 5. The first parameter of that distribution (as signified by the "1" between the two distribution names) follows a uniform distribution from 10 to 20. As another example, InfMix(Normal(0,5),2,Uniform(10,20)) defines a random variable that comes from a normal distribution with mean zero and standard deviation varying uniformly from 10 to 20. In general, the format of this distribution is: InfMix(ParentDist(Parms),MixParm,ParmDist(Parms)) where ParentDist is a distribution, MixParm is an integer indicating whether the first, second, ..., parameter of the ParentDist varies randomly, and ParmDist is the distribution of that parameter.

InfMix may be used recursively. For example, InfMix(InfMix(Normal(0,5),1,Uniform(0,2)),2,Uniform(4,6)) defines a normal distribution in which the mean is uniform(0,2) and the standard deviation is uniform(4,6).

Limitations: (1) At present, computations of the upper and lower bounds of InfMix distributions assume that the largest and smallest values of the random variable are obtained when the underlying ParmDist is at its two extremes. (2) Extreme caution is needed with these distributions because problems often arise in numerical integration. I have found it helpful to increase the IntegralMinSteps to 10, which was enough in most of the cases I've looked at, but you may need to adjust this up (for precision) or down (for speed) in your cases.

Truncated(RV,Min,Max)

A truncated distribution is a conditional distribution, conditioning on the random variable RV falling within the interval from Min to Max. For example, Truncated(Normal(0,1),-1,1) defines a random variable that is always between -1 and 1, and which within that interval has relative probabilities defined by the PDF of the standard normal. In general, the format of this distribution is: Truncated(BasisDistribution(Parms),Min,Max)

It is sometimes convenient to specify the truncation boundaries in terms the probabilities you want to cut off rather than the scores themselves. For example, you might want to look at the middle 90% of a normal distribution but might not immediately know which scores cut off the top and bottom 5%. For this reason, there is a variant of the command that takes probabilities instead of values for min and max, like this: TruncatedP(BasisDistribution(Parms),0.05,0.95) With TruncatedP, DistFit will use its InverseCDF function to find the score values that correspond to the cumulative probabilities that you specify, and then truncate at those score values.

Bounded(RV,Min,Max)

I do not know if this is a standard type of distribution or not, and would appreciate any comments on it from those in the know. A bounded distribution is similar to a truncated distribution in that the random variable must fall within the range of Min to Max. The difference is that all values less than Min are converted to Min, and all values less than Max are converted to Max. Thus, there are discrete masses of probability at Min and Max, and the probability density function between Min and Max is not conditionalized.

For example, consider the distribution Bounded(Normal(0,1),-1,1). This is really a mixture of these three distributions:

Distribution	Mixture Probability
Constant(-1)	0.1587
Truncated(Normal(0,1),-1,1)	0.6826
Constant(1)	0.1587

Note that 0.1587 is the probability that a normal(0,1) score is less than -1, and also the probability that it is greater than 1. Bounding the distribution thus means taking all of the probability density higher than the upper value and massing it at that value.

As with the truncated distribution, there is a form of the Bounded distribution based on probabilities, as in: BoundedP(Normal(0,1),0.1,0.9) .

Order(k,RV1,RV2,RV2,...,RVn)

The distribution of this order statistic is the distribution of the k'th largest observation in a sample of n independent observations from the n indicated random variables. For example, Order(2,Normal(0,1),Uniform(0,1),Exponential(1)) defines a random variable that is the median (2nd largest) in a sample containing one score from the standard normal, one from the uniform from 0-1, and one from the exponential with rate 1. In general, the format of this distribution is: Order(k,BasisDist1(Parms),...,BasisDistN(Parms)) In the special case where the basis distributions are all identical, it is more convenient to use the OrderIID distribution, described next.

OrderIID(k,n,RV)

This is the special case of the order distribution in which the basis distributions are identical as well as independent. In general, the format of this distribution is: OrderIID(k,N,BasisDist(Parms)) It is only necessary to specify the basis distribution once, since all are identical; instead, you have to specify how many there are (N).

OrdExp(i,n,l)

This is the special case of OrderIID in which the basis distribution is an exponential with rate l, and you want the i'th order statistic in a sample of n (1 £ i £ n). For this case there are nice fast closed forms for the mean and variance that were given to me by Rolf Ulrich.

OrdBinary(i,n1,RV1,n2,RV2)

This is an order distribution with two types of underlying RVs. For example, OrdBinary(2,5,Normal(0,1),7,Uniform(0,1)) specifies the distribution of the second order statistic in samples of 12 made up of five standard normals and seven standard uniforms.

MinBound(RV1,RV2)

Consider two arbitrary random variables X and Y, which may or may not be independent, and let Z º min(X,Y). The CDFs of these three random variables must obey the inequality

Fz(t) £ Fx(t) + Fy(t)     for all  t

because

Fz(t) = Fx(t) + Fy(t) - Pr (X £ t & Y £ t)

Thus, for any two basis RVs X and Y, we can construct the random variable Z which is a lower bound on the distribution of min(X,Y):

Fz(t) = ì ï í ï î Fx(t) + Fy(t) if  Fx(t) + Fy(t) < 1 1 if  Fx(t) + Fy(t) ³ 1

MinBound implements this lower bound distribution for any two arbitrary random variables X and Y.

8.5  Approximation Distributions

ApprPolygon(RV,NBins)

This is a continuous approximation that can be used only for a continuous basis distribution, RV. In brief, the PDF of the approximation distribution is a set of NBins straight lines, matched to the height of the basis distribution's PDF at the bin's edges (further detail is given below). For example, ApprPolygon(Convolution(Normal(0,1),Beta(2,2)),201) approximates the specified convolution with a set of 201 straight lines.

ApprFreqPolygon(RV,NBins)

This is a continuous approximation that can be used only for a continuous basis distribution, RV. In brief, the PDF of the approximation distribution is a set of NBins straight lines, matched to the height of the basis distribution's PDF at the bin's centers. For example, ApprFreqPolygon(Convolution(Normal(0,1),Beta(2,2)),201) approximates the specified convolution with a set of 201 straight lines.

This is my preferred type of approximation. It is generally quite accurate, and it is often much faster than the other approximations.

Details of construction.

Step 1:

The first line starts at X₁=minimum (of the basis distribution) with height PDF at that point and goes to X₂=minimum+W/2 with height PDF=Basis.PDF(X₂). The second line continues from the end of the first line to the point with X₃=minimum+1.5*W and height PDF=Basis.PDF(X₃). And so on, with the final line segment ending at the maximum of the basis distribution and PDF at the maximum.

Step 2:

The PDF just constructed is integrated, and the heights are scaled up or down appropriately so that the total area is 1.00.

ApprHistogram(RV,NBins)

This is a continuous approximation that can be used for either a discrete or a continuous basis distribution, RV. In brief, the PDF is of the approximation distribution is a set of NBins flat lines, as if the basis distribution were uniform within each bin (like in a histogram). For example, ApprHistogram(Convolution(Normal(0,1),Beta(2,2)),201) approximates the specified convolution with a set of 201 bins with equal probability within each bin.

This approximation is more general than ApprPolygon, because it can be used with discrete distributions, and it is less sensitive to abruptly-changing PDFs. But it is usually slower to construct initially, and it is often much slower to do any computations with. The PDF has discontinuities at the bin boundaries (unlike ApprPolygon), and these make numerical integrations converge more slowly.

Details of construction. The CDF of the basis distribution is computed at the top and bottom of each bin, and from these the bin probability is computed. Then, the height of the uniform approximation PDF within that bin is adjusted to give this bin probability.

ApprBinCen(RV,NBins)

This is a discrete approximation, and it can be used to approximate either a discrete or a continuous basis distribution, RV. In brief, the approximation assumes that all of the probability mass is concentrated in a single point at the center point of each bin; moreover, any value in the bin is treated as if it were that center point.

Details of construction. The CDF of the basis distribution is computed at the top and bottom of each bin, and from these the bin probability is computed. All of this probability mass is assigned to the value at the center of the bin. For purposes of PDF and CDF computations, all values within a bin are treated as equivalent to the center.

8.6  Distributions Arising in Connection with Signal Detection Theory

ZfromP(SampleSize,TrueP,Adjust)

This is the discrete distribution of Z, which is derived from the binomial distribution as follows:

1. For any sample from a Binomial(N,P), convert the number of successes k to the probability of success, p º k/N. If p=0, set p = Adjust/N; if p=1, set p=1-Adjust/N. "Adjust" is a parameter between 0 and 1, specified by the user, to indicate how the extreme data values should be treated.
2. Find Z such that p = Pr(z £ Z), where z is a random variable having the standard normal distribution.

APrime(NSignalTrials,PrHit,NNoiseTrials,PrFalseAlarm)

This is the distribution of the sample A¢ computed from an experiment with NSignalTrials signal trials each having the specified true probability of a hit, and NNoiseTrials noise trials each having the specified true probability of a false alarm. Specifically, A¢ is the distribution-free estimate of the area under the ROC curve computed using Equations 2 and 9 of Aaronson and Watts (1987).

APrimeSym(NTrials,PC)

This is a shortcut for the previous distribution that can be used when there are equal numbers of signal and noise trials and when the probability of a correct response (hit or correct rejection) is the same for both signal and noise trials.

YNdPrime(NSignalTrials,PrHit,NNoiseTrials,PrFalseAlarm,Adjust)

This is the distribution of the sample [^d]¢ computed from an experiment with NSignalTrials signal trials each having the specified true probability of a hit, NNoiseTrials noise trials each having the specified true probability of a false alarm, and using the Adjust factor (between 0 and 1) to correct cases with 0% or 100% hits or false alarms (e.g., replace 0 hits with Adjust hits, and replace NSignalTrials hits with [NSignalTrials - Adjust] hits). Programming note: If PDFs are requested, this distribution is implemented using the List (smaller samples) and AppApprCen (larger samples) random variables.

YNdPrimeSym(NTrials,TrueDP,Adjust)

This is the special case of YNdPrime in which NSignalTrials = NNoiseTrials and Pr(Hit) = 1 - Pr(FA). Note that the second parameter is the true d¢ rather than the hit probability.

9  Release History

10  Related Programs

Cupid

Interactive computations (pdf, cdf, moments, etc) with probability distributions. Can be used (for example) as an on-line table for distributions.

RandGen

Generates random values from a specified probability distribution.

DistFit

Estimates best-fitting parameters of a given probability distribution for a given data set.

MixTest

Computes a likelihood ratio test to see whether the difference between two conditions (say "experimental" versus "control") is a "uniform effect" or a "mixture effect". With a uniform effect, all of the scores in the experimental condition are increased relative to what they would have been in the control condition. With a mixture effect, however, only some of the scores in the experimental condition are affected; the rest of the scores in this condition are the same as they would have been without the manipulation (i.e., the same as they would have been in the control condition).

Pmetric

Estimates the parameters of a probability distribution from a data set relating the proportion of a certain (binary) response to a physical quantity. This type of analysis is often called "probit" analysis, and it is used (for example) in bioassay (analysis of dose/response curves) and psychophysics (analysis of psychometric functions).

11  Author Contact Address

12  Acknowledgements

13  References

Aaronson, D., & Watts, B. (1987). Extensions of Grier's computational formulas for A' and B" to below-chance performance. Psychological Bulletin, 102, 439-442.

D'Agostino, R. B. (1970). Simple compact portable test of normality: Geary's test revisited. Psychological Bulletin, 74, 138-140.

Dallal, G. E., & Wilkinson, L. (1986). An analytic approximation to the distribution of Lilliefors's test statistic for normality. The American Statistician, 40, 294-296.

Devroye, L. (1986). Non-uniform random variate generation. Berlin: Springer-Verlag.

Evans, M., Hastings, N., & Peacock, B. (1993). Statistical distributions. (2nd ed.) New York: Wiley.

Finney, D. J. (1978). Statistical method in biological assay. London: Griffin.

Gescheider, G. A. (1997). Psychophysics: The fundamentals. [3rd ed.]. Hillsdale, NJ: Erlbaum.

Johnson, N. L., & Kotz, S. (1970). Continuous univariate distributions. New York: Houghton Mifflin.

Johnson, N. L., Kotz, S., & Balakrishnan, N. (1994). Continuous univariate distributions. New York: Wiley.

Luce, R. D. (1986). Response times: Their role in inferring elementary mental organization. Oxford: Oxford University Press.

Quick, R. F. (1974). A vector magnitude model of contrast detection. Kybernetik, 16, 65-67.

Rosenbrock, H. H. (1960). An automatic method for finding the greatest or least value of a function. Computer Journal, 3, 175-184.

Schwarz, W. (2001). The ex-Wald distribution as a descriptive model of response times. Behavior Research Methods, Instruments, & Computers, 33, 457-469.

Schwarz, W. (2002). On the convolution of inverse Gaussian and exponential random variables. Communications in Statistics: Theory and Methods, 31, 2113-2121.

Sternberg, S., & Knoll, R. L. (1973). The perception of temporal order: Fundamental issues and a general model. In S. Kornblum (Ed.), Attention and performance IV (pp. 629-685). New York: Academic Press.

Strasburger, H. (2001). Converting between measures of slope of the psychometric function. Perception & Psychophysics, 63, 1348-1355.

14  Software License

14.1  Preamble

14.2  Terms of License

a) You must cause the modified files to carry prominent notices stating that you changed the files and the date of any change.

b) You must cause any work that you distribute or publish, that in whole or in part contains or is derived from the Program or any part thereof, to be licensed as a whole at no charge to all third parties under the terms of this License.

c) If the modified program normally reads commands interactively when run, you must cause it, when started running for such interactive use in the most ordinary way, to print or display an announcement including an appropriate copyright notice and a notice that there is no warranty (or else, saying that you provide a warranty) and that users may redistribute the program under these conditions, and telling the user how to view a copy of this License. (Exception: if the Program itself is interactive but does not normally print such an announcement, your work based on the Program is not required to print an announcement.)

a) Accompany it with the complete corresponding machine-readable source code, which must be distributed under the terms of Sections 1 and 2 above on a medium customarily used for software interchange; or,

b) Accompany it with a written offer, valid for at least three years, to give any third party, for a charge no more than your cost of physically performing source distribution, a complete machine-readable copy of the corresponding source code, to be distributed under the terms of Sections 1 and 2 above on a medium customarily used for software interchange; or,

c) Accompany it with the information you received as to the offer to distribute corresponding source code. (This alternative is allowed only for noncommercial distribution and only if you received the program in object code or executable form with such an offer, in accord with Subsection b above.)

Footnotes: