PaulTest4.java

package org.freehep.math.minuit.example.sim;

import java.util.ArrayList;
import java.util.List;
import java.util.Scanner;
import org.freehep.math.minuit.FCNBase;
import org.freehep.math.minuit.FunctionMinimum;
import org.freehep.math.minuit.MnMigrad;
import org.freehep.math.minuit.MnUserParameters;
import org.freehep.math.minuit.MnSimplex;

/**
 *
 * @version $Id: PaulTest4.java 8584 2006-08-10 23:06:37Z duns $
 */
public class PaulTest4
{
   public static void main(String[] args)
   {
      List<Double> positions = new ArrayList<Double>();
      List<Double> measurements = new ArrayList<Double>();
      List<Double> var = new ArrayList<Double>();
      
      Scanner in = new Scanner(PaulTest.class.getResourceAsStream("paul4.txt"));
      
      // read input data
      {
         while(in.hasNextDouble())
         {
            double x = in.nextDouble();
            double y = in.nextDouble();
            double err = in.nextDouble();
            
            positions.add(x);
            measurements.add(y);
            var.add(err*err);
         }
         System.out.printf("size= %d\n",var.size());
      }
      double[] m = new double[var.size()];
      double[] p = new double[var.size()];
      double[] v = new double[var.size()];
      for (int i=0; i<var.size(); i++)
      {
         m[i] = measurements.get(i);
         p[i] = positions.get(i);
         v[i] = var.get(i);
      }
      {
         // create Chi2 FCN function
         System.out.println(">>> test Chi2");
         PowerLawChi2FCN theFCN = new PowerLawChi2FCN(m, p, v);
         
         MnUserParameters upar = new MnUserParameters();
         upar.add("p0", -2.3, 0.2);
         upar.add("p1", 1100., 10.);
         
         MnMigrad migrad = new MnMigrad(theFCN, upar);
         migrad.setErrorDef(0.5);
         System.out.println("start migrad ");
         FunctionMinimum min = migrad.minimize();
         if(!min.isValid())
         {
            //try with higher strategy
            System.out.println("FM is invalid, try with strategy = 2.");
            migrad = new MnMigrad(theFCN, upar, 2);
            min = migrad.minimize();
         }
         System.out.println("minimum: "+min);
      }
      {
         System.out.println(">>> test log LikeliHood");
         // create LogLikelihood FCN function
         PowerLawLogLikeFCN theFCN = new PowerLawLogLikeFCN(m, p);
         
         MnUserParameters upar = new MnUserParameters();
         upar.add("p0", -2.1, 0.2);
         upar.add("p1", 1000., 10.);
         
         MnMigrad migrad = new MnMigrad(theFCN, upar);
         System.out.println("start migrad ");
         FunctionMinimum min = migrad.minimize();
         if(!min.isValid())
         {
            //try with higher strategy
            System.out.println("FM is invalid, try with strategy = 2.");
            migrad = new MnMigrad(theFCN, upar, 2);
            min = migrad.minimize();
         }
         System.out.println("minimum: "+min);
      }
      {
         System.out.println(">>> test Simplex");
         PowerLawChi2FCN chi2 = new PowerLawChi2FCN(m, p, v);
         PowerLawLogLikeFCN mlh = new PowerLawLogLikeFCN(m, p);
         
         MnUserParameters upar;
         double[] par = {-2.3, 1100.};
         double[] err = { 1., 1.};
         
         MnSimplex simplex = new MnSimplex(chi2, par, err);
         
         System.out.println("start simplex");
         FunctionMinimum min = simplex.minimize();
         System.out.println("minimum: "+min);
         
         MnSimplex simplex2 = new MnSimplex(mlh, par, err);   
         simplex2.setErrorDef(0.5);
         FunctionMinimum min2 = simplex2.minimize();
         System.out.println("minimum: "+min2);
      }
   }
   static class PowerLawFunc
   {
      
      PowerLawFunc(double p0, double p1)
      {
         theP0 = p0;
         theP1 = p1;
      }
      
      double valueOf(double x)
      {
         return p1()*Math.exp(Math.log(x)*p0());
      }
      
      double p0()
      {
         return theP0;
      }
      double p1()
      {
         return theP1;
      }
      
      private double theP0;
      private double theP1;
   };
   static class PowerLawChi2FCN implements FCNBase
   {
      
      PowerLawChi2FCN(double[] meas, double[] pos, double[] mvar)
      {
         theMeasurements = meas;
         thePositions = pos;
         theMVariances = mvar;
      }
      
      public double valueOf(double[] par)
      {
         assert(par.length == 2);
         PowerLawFunc pl = new PowerLawFunc(par[0], par[1]);
         double chi2 = 0.;
         
         for(int n = 0; n < theMeasurements.length; n++)
         {
            chi2 += ((pl.valueOf(thePositions[n]) - theMeasurements[n])*(pl.valueOf(thePositions[n]) - theMeasurements[n])/theMVariances[n]);
         }
         
         return chi2;
      }
      
      private double[] theMeasurements;
      private double[] thePositions;
      private double[] theMVariances;
   };
   static class PowerLawLogLikeFCN implements FCNBase
   {
      
      PowerLawLogLikeFCN(double[] meas, double[] pos)
      {
         theMeasurements = meas;
         thePositions = pos;
      }
      
      
      public double valueOf(double[] par)
      {
         assert(par.length == 2);
         PowerLawFunc pl = new PowerLawFunc(par[0], par[1]);
         double logsum = 0.;
         
         for(int n = 0; n < theMeasurements.length; n++)
         {
            double k = theMeasurements[n];
            double mu = pl.valueOf(thePositions[n]);
            logsum += (k*Math.log(mu) - mu);
         }
         
         return -logsum;
      }
      
      private double[] theMeasurements;
      private double[] thePositions;
   };
}