ROOT logo
#include "TRandom.h"
#include "TGraphErrors.h"
#include "TF1.h"
#include "TCanvas.h"
#include "TLegend.h"

void fitLinearRobust()
{
//This tutorial shows how the least trimmed squares regression,
//included in the TLinearFitter class, can be used for fitting
//in cases when the data contains outliers.
//Here the fitting is done via the TGraph::Fit function with option "rob":
//If you want to use the linear fitter directly for computing
//the robust fitting coefficients, just use the TLinearFitter::EvalRobust
//function instead of TLinearFitter::Eval
//Author: Anna Kreshuk
   
   //First generate a dataset, where 20% of points are spoiled by large
   //errors
   Int_t npoints = 250;
   Int_t fraction = Int_t(0.8*npoints);
   Double_t *x = new Double_t[npoints];
   Double_t *y = new Double_t[npoints];
   Double_t *e = new Double_t[npoints];
   TRandom r;
   Int_t i;
   for (i=0; i<fraction; i++){
      //the good part of the sample
      x[i]=r.Uniform(-2, 2);
      e[i]=1;
      y[i]=1 + 2*x[i] + 3*x[i]*x[i] + 4*x[i]*x[i]*x[i] + e[i]*r.Gaus();
   }
   for (i=fraction; i<npoints; i++){
      //the bad part of the sample
      x[i]=r.Uniform(-1, 1);
      e[i]=1;
      y[i] = 1 + 2*x[i] + 3*x[i]*x[i] + 4*x[i]*x[i]*x[i] + r.Landau(10, 5);
   }

   TGraphErrors *grr = new TGraphErrors(npoints, x, y, 0, e);
   grr->SetMinimum(-30);
   grr->SetMaximum(80);
   TF1 *ffit1 = new TF1("ffit1", "pol3", -5, 5);
   TF1 *ffit2 = new TF1("ffit2", "pol3", -5, 5);
   ffit1->SetLineColor(kBlue);
   ffit2->SetLineColor(kRed);
   TCanvas *myc = new TCanvas("myc", "Linear and robust linear fitting");
   myc->SetFillColor(42);
   myc->SetGrid();
   grr->Draw("ap");
   //first, let's try to see the result sof ordinary least-squares fit:
   printf("Ordinary least squares:\n");
   grr->Fit(ffit1);
   //the fitted function doesn't really follow the pattern of the data
   //and the coefficients are far from the real ones

   printf("Resistant Least trimmed squares fit:\n");
   //Now let's try the resistant regression
   //The option "rob=0.75" means that we want to use robust fitting and
   //we know that at least 75% of data is good points (at least 50% of points
   //should be good to use this algorithm). If you don't specify any number
   //and just use "rob" for the option, default value of (npoints+nparameters+1)/2
   //will be taken
   grr->Fit(ffit2, "+rob=0.75");
   //
   TLegend *leg = new TLegend(0.6, 0.8, 0.89, 0.89);
   leg->AddEntry(ffit1, "Ordinary least squares", "l");
   leg->AddEntry(ffit2, "LTS regression", "l");
   leg->SetFillColor(42);
   leg->Draw();

   delete [] x;
   delete [] y;
   delete [] e;

}
 fitLinearRobust.C:1
 fitLinearRobust.C:2
 fitLinearRobust.C:3
 fitLinearRobust.C:4
 fitLinearRobust.C:5
 fitLinearRobust.C:6
 fitLinearRobust.C:7
 fitLinearRobust.C:8
 fitLinearRobust.C:9
 fitLinearRobust.C:10
 fitLinearRobust.C:11
 fitLinearRobust.C:12
 fitLinearRobust.C:13
 fitLinearRobust.C:14
 fitLinearRobust.C:15
 fitLinearRobust.C:16
 fitLinearRobust.C:17
 fitLinearRobust.C:18
 fitLinearRobust.C:19
 fitLinearRobust.C:20
 fitLinearRobust.C:21
 fitLinearRobust.C:22
 fitLinearRobust.C:23
 fitLinearRobust.C:24
 fitLinearRobust.C:25
 fitLinearRobust.C:26
 fitLinearRobust.C:27
 fitLinearRobust.C:28
 fitLinearRobust.C:29
 fitLinearRobust.C:30
 fitLinearRobust.C:31
 fitLinearRobust.C:32
 fitLinearRobust.C:33
 fitLinearRobust.C:34
 fitLinearRobust.C:35
 fitLinearRobust.C:36
 fitLinearRobust.C:37
 fitLinearRobust.C:38
 fitLinearRobust.C:39
 fitLinearRobust.C:40
 fitLinearRobust.C:41
 fitLinearRobust.C:42
 fitLinearRobust.C:43
 fitLinearRobust.C:44
 fitLinearRobust.C:45
 fitLinearRobust.C:46
 fitLinearRobust.C:47
 fitLinearRobust.C:48
 fitLinearRobust.C:49
 fitLinearRobust.C:50
 fitLinearRobust.C:51
 fitLinearRobust.C:52
 fitLinearRobust.C:53
 fitLinearRobust.C:54
 fitLinearRobust.C:55
 fitLinearRobust.C:56
 fitLinearRobust.C:57
 fitLinearRobust.C:58
 fitLinearRobust.C:59
 fitLinearRobust.C:60
 fitLinearRobust.C:61
 fitLinearRobust.C:62
 fitLinearRobust.C:63
 fitLinearRobust.C:64
 fitLinearRobust.C:65
 fitLinearRobust.C:66
 fitLinearRobust.C:67
 fitLinearRobust.C:68
 fitLinearRobust.C:69
 fitLinearRobust.C:70
 fitLinearRobust.C:71
 fitLinearRobust.C:72
 fitLinearRobust.C:73
 fitLinearRobust.C:74
 fitLinearRobust.C:75
 fitLinearRobust.C:76
 fitLinearRobust.C:77