Defining Fitters for Java Analysis Studio

Right now, you'll have to look at how the least squares fitter is done and base your own on it.

package jas.hist.test;
import jas.hist.*;
import java.util.*;

public class LeastSquaresFit extends Fitter
{
    
    /* 
     * This algorithm based on CURFIT program from "Data Reduction and Error Analysis
     * for the Physical Sciences" by Philip R Bevington, Page 237. It performs a least
     * squared fit using the Marquardt algorithm.
     */

    public void fit(Fittable1DFunction func, double[] x, double[] y, double[] sigmaY) throws FitFailed
    {
        int npts = x.length;
        if (npts != y.length || npts != sigmaY.length) 
            throw new IllegalArgumentException("Invalid arguments to fit");

        double[] a = (double[]) func.getParameterValues().clone();
        int nterms = a.length;

        m_sigmaA = new double[nterms];

        int nFree = npts - nterms;
        if (nFree <= 0) throw new FitFailed("Insufficient degrees of Freedom");
        
        // Calculate weights

        double[] weight = new double[npts];
        for (int i=0; i<npts; i++) 
            if (sigmaY[i] == 0) weight[i] = 0; // ?????
            else weight[i] = 1/sigmaY[i]/sigmaY[i];

        double[] beta = new double[nterms];
        double[] b = new double[nterms];
        double[][] alpha = Matrix.create(nterms);
        double[][] array = Matrix.create(nterms);

        double lambda = 0.001;

        try
        {
            for (int iter = 0 ;; iter++)
            {
                if (iter > 10000) 
                    throw new FitFailed("Fit exceeded maximum number of iterations");         
            
                double chisqr = calculateChiSquared(func,x,y,weight,a,nFree);;

                // Evaluate Alpha and Beta matrices

                for (int j=0; j<nterms; j++)
                {
                    beta[j] = 0;
                    for (int k=0; k<=j; k++)
                    {
                        alpha[j][k] = 0;
                    }
                }

                for (int i=0; i<npts; i++)
                {
                    double[] deriv = func.getDerivatives(x[i],a);
                    for (int j=0; j<nterms; j++)
                    {
                        beta[j] += weight[i] * (y[i] - func.valueAt(x[i],a))*deriv[j];
                        for (int k=0; k<=j; k++)
                        {
                            alpha[j][k] += weight[i]*deriv[j]*deriv[k];
                        }
                    }
                }
                // Diagonalize matrix

                for (int j=0; j<nterms; j++)
                {
                    for (int k=0; k<j; k++)
                    {
                        alpha[k][j] = alpha[j][k];
                    }
                }

                // Evaluate chi squared at starting point

                double chisq1 = calculateChiSquared(func, x, y, weight, a, nFree);

                // Invert modified curvature matrix to find new parameters

                for ( ; ; lambda *= 10)
                {
                    for (int j=0; j<nterms; j++)
                    {
                        for (int k=0; k<nterms; k++)
                        {
                            array[j][k] = alpha[j][k] / Math.sqrt(alpha[j][j] * alpha[k][k]);
                        }
                        array[j][j] = 1 + lambda;
                    }
                    
                    try
                    {
                        double det = Matrix.invert(array);
                    }
                    catch (IndeterminateMatrix e)
                    {
                        throw new FitFailed("Matrix became indeterminate during fit");
                    }
                    
                    for (int j=0; j<nterms; j++)
                    {
                        b[j] = a[j];
                        for (int k=0; k<nterms; k++)
                        {
                            b[j] += beta[k]*array[j][k]/Math.sqrt(alpha[j][j]*alpha[k][k]);
                        }
                    }
                    
                    // If chi squared increased, increase lambda and try again
                    
                    chisqr = calculateChiSquared(func,x,y,weight,b,nFree);
                    if (chisqr <= chisq1) break;
                }
                
                // Evaluate parameters and uncertainties

                for (int j=0; j<nterms; j++)
                {
                    a[j] = b[j];
                    m_sigmaA[j] = Math.sqrt(array[j][j] / alpha[j][j]);
                }
                lambda /= 10;

                double delta = Math.abs(chisqr - m_chiSquared)/m_chiSquared;
                m_chiSquared = chisqr;
                if (delta < 0.001) break;

                setPercentComplete(10*iter);
            }
            try
            {
                func.setFit(this,a);
            }
            catch (InvalidFunctionParameter e)
            {
                throw new FitFailed("Unexpected exception: "+e);
            }
        }
        catch (FunctionValueUndefined e)
        {
            throw new FitFailed("Fit entered region where function value was undefined");
        }
    }
    private double calculateChiSquared(Fittable1DFunction func, double[] x, double[] y, double[] weight, double[] a, int nFree) throws FunctionValueUndefined
    {
        double result = 0;
        int npts = x.length;

        for (int i=0; i<npts; i++)
        {
            double delta = y[i] - func.valueAt(x[i],a);
            result += weight[i] * delta*delta;
        }
        return result/nFree;
    }
    public double getChiSquared()
    {
        return m_chiSquared;
    }
    public double[] getParameterSigmas()
    {
        return m_sigmaA;
    }
    private double m_chiSquared;
    private double[] m_sigmaA;
}

Page maintained by Jonas Gifford. Last updated 01/14/04.