ML_DISTFIT
Name
ML_DISTFIT
Purpose
Performs maximum likelihood fitting of a distribution.
Category
Math
Calling Sequence
ML_DISTFIT, X, Parm, Function_Name, ConfRegion
Inputs
X: Array of input data values. This is passed straight to the
user-supplied function, so complicated data structures that
encompass multi-dimensional information for each data
point can be used.
Parm: Variable containing initial guesses for parameters on input
and best fit values on output.
Function_Name: Name of user-supplied function defining the distribution.
The function must accept 2 arguments, X and Parm, and
return a vector containing the likelihood values for
each data point in X for the point in parameter space
given by Parm. The likelihood must be normalized so
that its total integral over all possible values of X
is a constant, regardless of Parm (it makes the most
sense to normalize this integral to unity, but that
is not strictly required).
Optional Outputs
ConfRegion: Lower and upper error estimates of each parameter,
marginalized over the other parameters.
I.e. ConfRegion[*,0] returns [low0,high0]
where low0 <= parm[0] <= high0
Keyword Parameters
FITA: Vector of which parameters should be fit (1 for each
parameter to be fit, 0 for each parameter to be held
constant).
THERE IS A BUG IN THE IMPLEMENTATION. DO NOT USE.
CONSTRAINT: Name of a user-supplied function that takes a parameter
vector as input and returns 1 if the point in parameter
space is permitted and 0 if it is not permitted.
LIKELIHOOD: Outputs an M-dimensional array with the likelihood
values over the range of parameter space probed. M is
the number of parameters that are fitted, which can be
less than the length of Parm if FITA is used.
LIKERANGE: 2xM dimensional array containing the bounds of the
LIKELIHOOD array.
Example
Fit the width and offset of a zero-centered Gaussian plus constant
distribution.
First, define the distribution function:
FUNCTION gauss_plus_const, X, Parm
; Parm[0]=constant offset, Parm[1]=width sigma
vmax = 2000.
normalization = Parm[1]*SQRT(!pi/2.)*ERF(vmax/(SQRT(2.)*Parm[1])) $
+ vmax*Parm[0]
distribution = EXP(-X^2/(2.*Parm[1]^2)) + Parm[0]
RETURN, distribution/normalization
END
Then generate some data that should adhere to this distribution,
with a width of 250 and a constant term containing 10% of the points.
IDL> data = [250*RANDOMN(seed, 900), 4000*(RANDOMU(seed, 100) - 0.5)]
And finally fit the distribution:
IDL> parm = [0., 100.]
IDL> ML_DISTFIT, data, parm, 'gauss_plus_const', parmconf
IDL> PRINT, parm
0.0625345 223.94577
IDL> PRINT, parmconf
0.057511609 0.0973332
207.087 243.841
Modification History
Written by: Jeremy Bailin. Thanks to the writers of MLEfit.pro, which
furnished the Hessian routines, Peder Norberg for useful
discussions, and Nicolas Petitclerc for additional testing.
27 Nov 2008 Release in JBIU.