Fitter

The AIDA IFitter interface provides the user the possibility to fit IFunctions to any AIDA data storage object. Binned fits can be perfomed on IHistograms, IProfiles and IDataPointSets, while unbinned fits can be performed on IClouds and ITuples. Simple fits can be perfomed directely on the data storage objects while the IFitData interface is to be used for a greater control over the data, in particular its ranges and the connection to the IFunction's variables. Through the IFitter it is also possible to change the underlying optimization engine as well as the fit method used.

Chi2 Fit to an IHistogram

A simple chi2 fit to an histogram is perfomed in the example below:

import hep.aida.*;
import java.util.Random;

public class Chi2FitToHistogram
{
   public static void main(String[] args)
   {
      // Create factories
      IAnalysisFactory analysisFactory = IAnalysisFactory.create();
      IHistogramFactory histogramFactory = analysisFactory.createHistogramFactory(analysisFactory.createTreeFactory().create());
      IPlotter plotter = analysisFactory.createPlotterFactory().create("Plot");
      IFitFactory fitFactory = analysisFactory.createFitFactory();
      
      // Create 1D histogram
      IHistogram1D h1d = histogramFactory.createHistogram1D("Gaussian Distribution",100,-5,5);
      
      // Fill 1D histogram with Gaussian
      Random r = new Random();
      for (int i=0; i<5000; i++)
         h1d.fill(r.nextGaussian());
      
      // Do Fit
      IFitter fitter = fitFactory.createFitter("chi2");
      IFitResult result = fitter.fit(h1d,"g");
      
      // Show results
      plotter.createRegions(1,1,0);
      plotter.region(0).plot(h1d);
      plotter.region(0).plot(result.fittedFunction());
      plotter.show();
   }
}

Simple Binned and Unbinned Fits

In the example below it is shown how to perfom a binned and an unbinned fit over the same data using different optimizers and fit methods.

In the JAIDA implementation there are currently two optimizers available: Minuit and Uncmin and it is possible to choose among the following fit methods: ChiSquared, CleverChiSquared, BinnedMaximumLikelihood, LeastSquares, UnbinnedMaximumLikelihood.

import hep.aida.*;
import java.util.Random;

public class SimpleFit
{
  public static void main(String[] argv) throws java.io.IOException
  {
    
    IAnalysisFactory  anFactory   = IAnalysisFactory.create();
    ITree             tree        = anFactory.createTreeFactory().create("storeForFittingExample");
    IHistogramFactory histFactory = anFactory.createHistogramFactory( tree );
    ITupleFactory     tuplFactory = anFactory.createTupleFactory( tree );
    IFunctionFactory  funcFactory = anFactory.createFunctionFactory( tree );
    IFitFactory       fitFactory  = anFactory.createFitFactory();
    IFitter           fitter      = fitFactory.createFitter("Chi2","minuit");
    
    IHistogram1D gaussHist  = histFactory.createHistogram1D("gaussHist","Gaussian Histogram",100,-5,5);
    ICloud1D     gaussCloud = histFactory.createCloud1D("gaussCloud","Gaussian Cloud");
    ITuple       tuple      = tuplFactory.create("tuple","Tuple Example","double gaussDistr");

    Random r = new Random();
    for (int i=0; i<10000; i++) {
      double x = r.nextGaussian();
      gaussHist.fill(x);
      gaussCloud.fill(x);
      tuple.fill(0,x);
      tuple.addRow();
    }

    // Chi2 fit with Minuit
    IFitResult minuitChi2Fit = fitter.fit(gaussHist,"g");
    
    // Least Squares fit with Minuit
    fitter.setFitMethod("LS");
    IFitResult minuitLeastSquaresFit = fitter.fit(gaussHist,"g");      

    // Binned Maximum Likelihood fit with Minuit
    fitter.setFitMethod("binnedMaximumLikelihood");
    IFitResult minuitBinnedMaxLikelihoodFit = fitter.fit(gaussHist,"g"); 

    // Unbinned Maximum Likelihood fit with Uncmin
    fitter.setEngine("uncmin");
    fitter.setFitMethod("uml");
    IFitResult uncminUMLFitToCloud = fitter.fit(gaussCloud,"g");

    String[] gaussColumn = {"gaussDistr"};
    IFitData fitData = fitFactory.createFitData();
    fitData.createConnection(tuple,gaussColumn);

    IFitResult uncminUMLFitToTuple = fitter.fit(fitData,"g");


    IHistogram1D gaussProj  = histFactory.createHistogram1D("gaussProj","Gaussian Histogram Projected from ITuple",100,-5,5);
    tuple.project( gaussProj, tuplFactory.createEvaluator("gaussDistr") );
    IPlotter plotter = anFactory.createPlotterFactory().create("Plot");
    plotter.createRegion(0,0,.66,1).plot(gaussHist);
    plotter.destroyRegions();
    plotter.createRegion(0,0,.66,1).plot(gaussHist);
    plotter.region(0).plot( minuitChi2Fit.fittedFunction() );
    plotter.region(0).plot( minuitLeastSquaresFit.fittedFunction() );
    plotter.region(0).plot( minuitBinnedMaxLikelihoodFit.fittedFunction() );
    
    gaussCloud.convert(100,gaussCloud.lowerEdge(),gaussCloud.upperEdge());
    plotter.createRegion(.66,0,.33,.5).plot(gaussCloud.histogram());
    
    IModelFunction cloudFunc = (IModelFunction)uncminUMLFitToCloud.fittedFunction();
    cloudFunc.normalize(false);
    double gaussCloudNorm = gaussCloud.entries()*( gaussCloud.upperEdge()-gaussCloud.lowerEdge() )/gaussCloud.histogram().axis().bins();
    double cloudFuncNorm = gaussCloudNorm/(Math.sqrt(2*Math.PI)*cloudFunc.parameter("sigma"));
    cloudFunc.setParameter("amplitude",cloudFuncNorm);
    plotter.region(1).plot( cloudFunc );

    plotter.createRegion(.66,.5,.33,.5).plot( gaussProj );


    IModelFunction tupleFunc = (IModelFunction)uncminUMLFitToTuple.fittedFunction();
    tupleFunc.normalize(false);
    
    double gaussProjNorm = gaussProj.entries()*( gaussProj.axis().upperEdge()-gaussProj.axis().lowerEdge() )/gaussProj.axis().bins();
    double tupleFuncNorm = gaussProjNorm/(Math.sqrt(2*Math.PI)*tupleFunc.parameter("sigma"));

    tupleFunc.setParameter("amplitude",tupleFuncNorm);
    plotter.region(2).plot( tupleFunc );
    plotter.show();
  }
}

Control fit parameters and set constraints

With the IFitter interface it is possible to have a more direct control over the fit: the parameters in the fit can be controlled with the IFitParameterSettings interface and it is also possible to set constraints among them. In the example below we create a scripted function an show how to control its parameters:

import hep.aida.*;
import java.util.Random;

public class ComplexFit
{
  public static void main(String[] argv) throws java.io.IOException
  {
    
    IAnalysisFactory  anFactory   = IAnalysisFactory.create();
    ITree             tree        = anFactory.createTreeFactory().create("storeForFittingExample");
    IHistogramFactory histFactory = anFactory.createHistogramFactory( tree );
    ITupleFactory     tuplFactory = anFactory.createTupleFactory( tree );
    IFunctionFactory  funcFactory = anFactory.createFunctionFactory( tree );
    IFitFactory       fitFactory  = anFactory.createFitFactory();
    IFitter           fitter      = fitFactory.createFitter("Chi2","minuit");
    
    IHistogram2D hist  = histFactory.createHistogram2D("hist","Test Histogram",100,-5,15,50,-5,5);

    Random r = new Random();

    for (int i=0; i<10000; i++) {
      double x = r.nextGaussian()+5;
      if ( r.nextDouble() > 0.8 ) x = 2.5*r.nextGaussian()+5;
      double y = r.nextGaussian();
      hist.fill(x,y);
    }
    
    IFitData fitData = fitFactory.createFitData();
    fitData.create2DConnection(hist);

    IFunction func = funcFactory.createFunctionFromScript("twoDdistr",2,"N*(a*exp( -(x[0]-mu0)*(x[0]-mu0)/(2*s0*s0) )+"
                       +"(1-a)*exp( -(x[0]-mu1)*(x[0]-mu1)/(2*s1*s1) ))*exp( -(x[1]-mu2)*"
                       +"(x[1]-mu2)/(2*s2*s2) )","N,a,mu0,s0,mu1,s1,mu2,s2","",null);


    double[] initialPars = { 1, 0.8, 5, 1, 5, 2, 0, 1};
    func.setParameters( initialPars );

    fitter.fitParameterSettings("mu2").setFixed(true);
    fitter.fitParameterSettings("a").setBounds(0.5,0.9);
    fitter.fitParameterSettings("a").setStepSize(0.001);
    fitter.fitParameterSettings("s1").setBounds(2,4);
    fitter.fitParameterSettings("s1").setStepSize(0.1);
    fitter.setConstraint("s0 = s2");
    fitter.setConstraint("mu0 = mu1");


    IFitResult fitResult = fitter.fit(fitData,func);

    System.out.println("Chi2 = "+fitResult.quality());
    
    double[] fPars     = fitResult.fittedParameters();
    double[] fParErrs  = fitResult.errors();
    String[] fParNames = fitResult.fittedParameterNames();

    for(int i=0; i< fitResult.fittedFunction().numberOfParameters(); i++ ) 
      System.out.println(fParNames[i]+" : "+fPars[i]+" +- "+fParErrs[i]);


    IPlotter plotter = anFactory.createPlotterFactory().create("Plot");
    plotter.destroyRegions();
    plotter.createRegion(0,0,.66,1).plot(hist);
    plotter.createRegion(.66,0,.33,.5).plot( histFactory.projectionX("projX",hist) );
    plotter.createRegion(.66,.5,.33,.5).plot( histFactory.projectionY("projY",hist) );
    plotter.show();

  }
}

Scans and contours

Scans and contours can be obtained as shown below:

import hep.aida.*;
import java.util.Random;

public class ScanAndContour
{       
  public static void main(String[] argv) throws java.io.IOException
  {        
    IAnalysisFactory  anFactory   = IAnalysisFactory.create();
    ITree             tree        = anFactory.createTreeFactory().create("storeForFittingExample");
    IHistogramFactory histFactory = anFactory.createHistogramFactory( tree );
    ITupleFactory     tuplFactory = anFactory.createTupleFactory( tree );
    IFunctionFactory  funcFactory = anFactory.createFunctionFactory( tree );
    IFitFactory       fitFactory  = anFactory.createFitFactory();
    IFitter           fitter      = fitFactory.createFitter("Chi2","minuit");

    IHistogram1D hist  = histFactory.createHistogram1D("hist","Test Histogram",100,-5,5);

    Random r = new Random();

    for (int i=0; i<10000; i++) {
      double x = r.nextGaussian();
      hist.fill(x);
    }

   IFitData fitData = fitFactory.createFitData();
   fitData.create1DConnection(hist);

   IFitResult fitResult = fitter.fit(fitData,"g");

   int meanIndex  = fitResult.fittedFunction().indexOfParameter("mean");
   double meanVal = fitResult.fittedParameters()[meanIndex];
   double meanErr = fitResult.errors()[meanIndex];
   IDataPointSet meanScan = fitter.createScan1D(fitData, fitResult.fittedFunction(),"mean",20, meanVal-3*meanErr, meanVal+3*meanErr);

   int sigmaIndex  = fitResult.fittedFunction().indexOfParameter("sigma");
   double sigmaVal = fitResult.fittedParameters()[sigmaIndex];
   double sigmaErr = fitResult.errors()[sigmaIndex];
   IDataPointSet oneSigmaContour = fitter.createContour(fitData, fitResult, "mean", "sigma", 10, 1);
   IDataPointSet twoSigmaContour = fitter.createContour(fitData, fitResult, "mean", "sigma", 10, 2);

   IPlotter plotter = anFactory.createPlotterFactory().create("Plot");
   plotter.destroyRegions();
   plotter.createRegion(0,0,.66,1).plot(hist);
   plotter.region(0).plot(fitResult.fittedFunction());
   plotter.createRegion(.66,0,.33,.5).plot( meanScan );
   plotter.createRegion(.66,.5,.33,.5).plot( twoSigmaContour );
   plotter.region(2).plot( oneSigmaContour );
   plotter.show();

  }
}

Table of Contents


$Id: fitter.jsp,v 1.1 2006/03/09 18:22:10 turri Exp $\