JMEstimator.hh

Go to the documentation of this file.

#ifndef __JFIT__JMESTIMATOR__
#define __JFIT__JMESTIMATOR__
 
#include <cmath>
 
 
/**
 * \file
 * Maximum likelihood estimator (M-estimators).
 * \author mdejong
 */
namespace JFIT {}
namespace JPP { using namespace JFIT; }
 
namespace JFIT {
 
  /**
   * Interface for maximum likelihood estimator (M-estimator).
   */
  struct JMEstimator {
    /**
     * Virtual destructor.
     */
    virtual ~JMEstimator()
    {}
 
 
    /**
     * Get maximum likelihood estimate.
     *
     * \param  z           deviation
     * \return             likelihood
     */
    virtual double getRho(const double z) const = 0;
 
 
    /**
     * Get derivative of maximum likelihood estimate.
     *
     * \param  z           deviation
     * \return             derivative
     */
    virtual double getPsi(const double z) const = 0;
  };
 
 
  /**
   * Null M-estimator.
   *
   * This is not an estimator at all, but an object that just returns whatever it is given.
   * It was introduced so that the user can directly access the likelihood calculated by JRegressor<JEnergy>.
   */
  struct JMEstimatorNull :
    public JMEstimator
  {
    virtual double getRho(const double z) const { return z; }
    virtual double getPsi(const double z) const { return 1.0; }
  };
 
 
  /**
   * Normal M-estimator.
   *
   * This estimator is based on a Gaussian PDF and produces the standard chi2.
   */
  struct JMEstimatorNormal :
    public JMEstimator
  {
    virtual double getRho(const double z) const { return 0.5*z*z; }
    virtual double getPsi(const double z) const { return z; }
  };
 
 
  /**
   * Lorentzian M-estimator.
   *
   * This estimator prouces a logarithmic dependence for large deviations.
   */
  struct JMEstimatorLorentzian :
    public JMEstimator
  {
    virtual double getRho(const double z) const { return log (1.0 + 0.5*z*z); }
    virtual double getPsi(const double z) const { return z / (1.0 + 0.5*z*z); }
  };
 
 
  /**
   * Linear M-estimator.
   *
   * This estimator produces a linear dependence for large deviations.
   */
  struct JMEstimatorLinear :
    public JMEstimator
  {
    virtual double getRho(const double z) const { return sqrt(1.0 + 0.5*z*z) - 1.0; }
    virtual double getPsi(const double z) const { return 0.5 * z / sqrt(1.0 + 0.5*z*z); }
  };
 
 
  /**
   * Tukey's biweight M-estimator.
   *
   * This estimator produces a redescending dependence for large deviations.
   */
  struct JMEstimatorTukey :
    public JMEstimator
  {
    /**
     * Constructor.
     *
     * \param  __k         standard deviation
     */
    JMEstimatorTukey(const double __k) :
      k(__k)
    {}
 
    virtual double getRho(const double z) const
    { 
      const double w = 0.5 * k*k / 3.0;
 
      if (fabs(z) < k) {
 
        const double u = z/k;
        const double v = 1.0 - u*u;
        
        return w * (1.0 - v*v*v);
      }
 
      return w;
    }
 
    virtual double getPsi(const double z) const 
    { 
      if (fabs(z) < k) {
 
        const double u = z/k;
        const double v = 1.0 - u*u;
        
        return z * v*v;
      }
 
      return 0.0;
    }
 
    double k;
  };
}
 
#endif