* Given a data series, a model function, and an initial guess for the parameter * values, class RobustFit's fit() method finds parameter values that * minimize the following metric: *
* Σi * ρ (yi − * f (xi, parameters)) *
* where f is the model function and ρ is one of these metric * functions: *
* In other words, the fit() method fits the model to the data * by adjusting the parameters to minimize the metric. *
* The metric function is the negative logarithm of the probability distribution * of the errors in the y values. The above metric functions correspond * to normal, two-sided exponential, and Cauchy error distributions. *
* The metric functions differ in how they treat outliers, i.e., data * points that deviate from the model. The normal metric function gives * increasing weights to points with increasing deviations. However, because of * the increasing weights, outlier points may skew the fit (hence, this is not * really a "robust" metric function). The exponential metric function gives * equal weights to all points, regardless of deviation. This reduces the * influence of outliers on the fit, yielding a more robust fit. With the * Cauchy metric function, the weights first increase, then decrease as the * deviations increase. This reduces the influence of outliers even further. *
* The fit() method uses class {@linkplain * MDMinimizationDownhillSimplex} to find the parameter values that minimize the * metric. The inputs to and outputs from the fit() method are stored * in fields of an instance of class RobustFit. *
* The fitWithDistribution() method uses the bootstrapping * technique to determine the distribution of the model parameters, which * depends on the error distribution of the data points. Bootstrapping performs * multiple iterations of the model fitting procedure. On each iteration, a * trial data set the same size as the original data set is created by sampling * the original data points with replacement, and model parameters for the trial * data set are computed. The fitWithDistribution() method outputs a * series of the parameter values found at each iteration; the confidence region * for the parameters; and the goodness-of-fit p-value. * * @author Alan Kaminsky * @version 22-Oct-2010 */ public class RobustFit { // Exported data members. /** * The model function. When model.f() is called, the x * argument is xi, the x value of a data * point; the p argument contains the model parameters; and the * return value is * f (xi, parameters). */ public final ParameterizedFunction model; /** * The number of parameters in the model, M. */ public final int M; /** * The metric function. By default, this is CAUCHY. It can instead * be set to NORMAL, EXPONENTIAL, or some other metric * function. */ public Function metric = CAUCHY; /** * The model parameters. On input to the fit() and * fitWithDistribution() methods, param contains the * initial guess for the model parameters. On output from the fit() * and fitWithDistribution() methods, param contains the * fitted parameter values. */ public final double[] param; /** * The data series. It contains the (x,y) data points to be fitted to * the model. It is specified as an argument of the fit() and * fitWithDistribution() methods. */ public XYSeries data; /** * The metric value. An output of the fit() and * fitWithDistribution() methods. It is set to the value of the * metric for the model with the fitted parameters stored in param. */ public double metricValue; /** * The model parameter distribution. An output of the * fitWithDistribution() method. paramSeries is a * T-element array, where T is the number of trials. Each * element of paramSeries is an M-element array giving the * fitted parameter values for the corresponding trial. */ public double[][] paramSeries; /** * The metric values for the model parameter distribution. An output of the * fitWithDistribution() method. metricSeries is a * T-element array, where T is the number of trials. Each * element of metricSeries gives the value of the metric for the * model with the parameters stored in the corresponding element of * paramSeries. */ public double[] metricSeries; /** * The lower bound of the confidence region for the model parameters. An * output of the fitWithDistribution() method. The confidence level * is specified as an argument of the fitWithDistribution() method; * for example, 0.90 specifies a 90% confidence level. The confidence region * is an M-dimensional rectangular hyperprism centered on the fitted * parameters stored in param, such that the given fraction of the * model parameter distribution stored in paramSeries falls within * the hyperprism. confidenceRegionLowerBound gives the lower bound * of each dimension of the confidence region hyperprism. */ public double[] confidenceRegionLowerBound; /** * The upper bound of the confidence region for the model parameters. An * output of the fitWithDistribution() method. * confidenceRegionUpperBound gives the upper bound of each * dimension of the confidence region hyperprism. */ public double[] confidenceRegionUpperBound; /** * The goodness-of-fit p-value. An output of the * fitWithDistribution() method. This gives the probability that a * metric value greater than or equal to metricValue would occur by * chance, even if the model with parameters params is correct. */ public double pValue; // Exported constants. /** * The normal metric function. */ public static final Function NORMAL = new Function() { public double f (double z) { return 0.5*z*z; } }; /** * The exponential metric function. */ public static final Function EXPONENTIAL = new Function() { public double f (double z) { return abs(z); } }; /** * The Cauchy metric function. */ public static final Function CAUCHY = new Function() { public double f (double z) { return log(1.0 + 0.5*z*z); } }; // Hidden data members. // For minimizing the metric. private MDMinimizationDownhillSimplex minimizer; // Hidden helper classes. /** * Class Metric computes the metric to be minimized. * * @author Alan Kaminsky * @version 06-Oct-2010 */ private class Metric implements MDFunction { public int argumentLength() { return M; // Number of model parameters } public double f (double[] p) // Model parameters { double sum = 0.0; int N = data.length(); for (int i = 0; i < N; ++ i) { sum += metric.f (data.y(i) - model.f (data.x(i), p)); } return sum; } } // Exported constructors. /** * Construct a new robust fitting object for the given model. The * model field is set to the corresponding argument. The M * field is set by calling the model function's parameterLength() * method. The param field is allocated with M elements; * initially, the elements are 0. * * @param model Model function. * * @exception NullPointerException * (unchecked exception) Thrown if model is null. */ public RobustFit (ParameterizedFunction model) { if (model == null) { throw new NullPointerException ("RobustFit(): model is null"); } this.model = model; this.M = model.parameterLength(); this.param = new double [M]; minimizer = new MDMinimizationDownhillSimplex (new Metric()); } // Exported operations. /** * Fit the given data series to the model. The data series is stored in the * data field. The model function was specified to the constructor, * and is also stored in the model field. On input to the * fit() method, param contains the initial guess for the * model parameters. On output from the fit() method, * param contains the fitted parameter values and * metricValue contains the value of the metric for the fitted * parameters. *
* The fit() method uses the downhill simplex technique to find the * model parameters that minimize the metric. This involves initializing the * simplex in an {@linkplain MDMinimizationDownhillSimplex} object. * The initializeSimplex() method is called to initialize the * simplex. * * @param data Data series. * * @exception TooManyIterationsException * (unchecked exception) Thrown if too many iterations occurred without * finding parameters that minimize the metric function. */ public void fit (XYSeries data) { this.data = data; initializeSimplex (minimizer); minimizer.minimize(); System.arraycopy (minimizer.x[0], 0, param, 0, M); metricValue = minimizer.f[0]; } /** * Fit the given data series to the model and compute the distribution of * the model parameters. The data series is stored in the data * field. The bootstrapping technique with T trials using the given * pseudorandom number generator is used to compute the distribution. The * given confidence level is used to compute the confidence region; for * example, 0.90 specifies a 90% confidence level. The model function was * specified to the constructor, and is also stored in the model * field. On input to the fitWithDistribution() method, * param contains the initial guess for the model parameters. On * output from the fit() method, param contains the fitted * parameter values, metricValue contains the value of the metric * for the fitted parameters, paramSeries contains the series of * fitted parameter values from all the trials, metricSeries * contains the metric values from all the trials, * confidenceRegionLowerBound and * confidenceRegionUpperBound contain the lower and upper bounds of * the confidence region hyperprism, and pValue contains the * goodness-of-fit. *
* The fitWithDistribution() method uses the downhill simplex
* technique to find the model parameters that minimize the metric. This
* involves initializing the simplex in an {@linkplain
* MDMinimizationDownhillSimplex} object. The initializeSimplex()
* method is called to initialize the simplex.
*
* @param data Data series.
* @param T Number of trials.
* @param prng Pseudorandom number generator.
* @param conf Confidence level, in the range 0.0 .. 1.0.
*
* @exception IllegalArgumentException
* (unchecked exception) Thrown if conf is out of bounds.
* @exception TooManyIterationsException
* (unchecked exception) Thrown if too many iterations occurred without
* finding parameters that minimize the metric function.
*/
public void fitWithDistribution
(XYSeries data,
int T,
Random prng,
double conf)
{
if (0.0 > conf || conf > 1.0)
{
throw new IllegalArgumentException
("RobustFit.fitWithDistribution(): conf = "+conf+" illegal");
}
int N = data.length();
// Allocate storage for outputs.
paramSeries = new double [T] [M];
metricSeries = new double [T];
confidenceRegionLowerBound = new double [M];
confidenceRegionUpperBound = new double [M];
// Compute fitted parameters for original data set.
this.data = data;
initializeSimplex (minimizer);
minimizer.minimize();
System.arraycopy (minimizer.x[0], 0, param, 0, M);
metricValue = minimizer.f[0];
// Perform bootstrap trials.
ListXYSeries trialData = new ListXYSeries();
this.data = trialData;
for (int trial = 0; trial < T; ++ trial)
{
// Create trial data set.
trialData.clear();
for (int i = 0; i < N; ++ i)
{
int j = prng.nextInt (N);
trialData.add (data.x(j), data.y(j));
}
// Perform robust fit on trial data set.
initializeSimplex (minimizer);
minimizer.minimize();
// Record fitted parameters and metric value.
System.arraycopy (minimizer.x[0], 0, paramSeries[trial], 0, M);
metricSeries[trial] = minimizer.f[0];
}
this.data = data;
// Compute confidence region: First, compute distance of each
// bootstrapped parameter set from the fitted parameter set.
final double[] distance = new double [T];
for (int trial = 0; trial < T; ++ trial)
{
double[] paramSeries_trial = paramSeries[trial];
double maxDistance = 0.0;
for (int i = 0; i < M; ++ i)
{
maxDistance =
max (maxDistance, abs (paramSeries_trial[i] - param[i]));
}
distance[trial] = maxDistance;
}
// Compute confidence region: Second, create an index array for visiting
// the distance, paramSeries, and metricSeries arrays in ascending order
// of distance.
Integer[] index = new Integer [T];
for (int i = 0; i < T; ++ i) index[i] = i;
Arrays.sort (index, new Comparator
* The default implementation of this method sets the first simplex point to
* param, sets the second simplex point to param except
* element 0 is set to perturb(param[0]), sets the third simplex
* point to param except element 1 is set to
* perturb(param[1]), and so on. perturb(x) = 1.01x
* if x ≠ 0; perturb(0) = 0.01. This method can be overridden to
* initialize the simplex differently.
*/
protected void initializeSimplex
(MDMinimizationDownhillSimplex minimizer)
{
System.arraycopy (param, 0, minimizer.x[0], 0, M);
for (int i = 1; i <= M; ++ i)
{
double[] x_i = minimizer.x[i];
System.arraycopy (param, 0, x_i, 0, M);
x_i[i-1] = perturb (x_i[i-1]);
}
}
private static double perturb
(double x)
{
return x == 0.0 ? 0.01 : 1.01*x;
}
}