Linear Regression

<< Graphics`Graphics`

First, let's create a data set which we would like to fit a straight line to

num = 20;
dataset  =  {{1, 7.18437, 1.13691}, {2, 7.41171, 1.41883}, {3, 
            6.46758, 1.46121}, {4, 9.30980, 1.36427}, {5, 7.54707, 
            1.18999}, {6, 12.6665, 1.75589}, {7, 12.63065, 1.586572}, {8, 
            11.73464, 1.009831}, {9, 14.28273, 1.277198}, {10, 15.95493, 
            1.617627}, {11, 16.86929, 1.846311}, {12, 15.17158, 
            1.479319}, {13, 19.68491, 1.709392}, {14, 18.75987, 
            1.060488}, {15, 20.21733, 1.248179}, {16, 22.45007, 
            1.696217}, {17, 21.67026, 1.058186}, {18, 24.78354, 
            1.940323}, {19, 25.03960, 1.471614}, {20, 25.04889, 
            1.930492}};

datasetplot  =  ErrorListPlot[dataset]

Create the sums that we need:

sum  =  
      Sum[1/Part[dataset, i, 3]^2, {i, 1, num}]
sumt  =  
      Sum[Part[dataset, i, 1]/Part[dataset, i, 3]^2, {i, 1, 
          num}]
sumy  =  
      Sum[Part[dataset, i, 2]/Part[dataset, i, 3]^2, {i, 1, 
          num}]
sumtt  =  
      Sum[Part[dataset, i, 1]^2/Part[dataset, i, 3]^2, {i, 1, 
          num}]
sumty  =  
      Sum[Part[dataset, i, 1]*Part[dataset, i, 2]/Part[dataset, i, 3]^2, {i, 
          1, num}]
delta  =  sum*sumtt  -  sumt^2

Create the linear regression model:

a  =  (sumtt*sumy - sumt*sumty)/
          delta;
b  =  (sum*sumty - sumt*sumy)/
          delta;
f[t_]  =  a  +  b*t
model  =  Plot[f[t], {t, 0, num}]

Show[model, datasetplot, Frame -> True, 
      FrameLabel -> {"t", "y"}]

So, what if we did not have the error bars? Let's see what we would have:

datasetreduced  =  {{1, 7.18437}, {2, 7.41171}, {3, 6.46758}, {4, 
            9.30980}, {5, 7.54707}, {6, 12.6665}, {7, 12.63065}, {8, 
            11.73464}, {9, 14.28273}, {10, 15.95493}, {11, 16.86929}, {12, 
            15.17158}, {13, 19.68491}, {14, 18.75987}, {15, 20.21733}, {16, 
            22.45007}, {17, 21.67026}, {18, 24.78354}, {19, 25.03960}, {20, 
            25.04889}};

datasetreducedplot  = ListPlot[datasetreduced]

sum  =  num
sumt  =  
      Sum[Part[datasetreduced, i, 1], {i, 1, num}]
sumy  =  
      Sum[Part[datasetreduced, i, 2], {i, 1, num}]
sumtt  =  
      Sum[Part[datasetreduced, i, 1]^2, {i, 1, 
          num}]
sumty  =  
      Sum[Part[datasetreduced, i, 1]*Part[datasetreduced, i, 2], {i, 1, 
          num}]
delta  =  sum*sumtt  -  sumt^2

Create the linear regression model:

a  =  (sumtt*sumy - sumt*sumty)/
          delta;
b  =  (sum*sumty - sumt*sumy)/
          delta;
f[t_]  =  a  +  b*t
model  =  Plot[f[t], {t, 0, num}]

compplot  =  
      Show[model, datasetreducedplot, Frame -> True, 
        FrameLabel -> {"t", "y"}]

Error Estimation

How good is our fit? We can calculate the correlation coefficient to find out.

ybar  =  
      Sum[Part[datasetreduced, i, 2], {i, 1, num}]/
        num
tbar  =  
      Sum[Part[datasetreduced, i, 1], {i, 1, num}]/
        num
Rsquare  =  
      Sum[(f[Part[datasetreduced, i, 1]]  -   ybar)^2, {i, 1, num}]/
        Sum[(Part[datasetreduced, i, 2] -   ybar)^2, {i, 1, 
            num}]

Since we calculate (R^2)to be about 0.96~1, we have confidence that the data actually represent a linear model, and the model represents the data well.

For the data set, we also can verify graphically that the model function goes through the
centroid of the data set.

center  =  
      Graphics[Disk[{tbar, ybar}, 0.4]]
Show[compplot, center, AspectRatio -> 1]

Extrapolation

Once we have the model, we can use it to extrapolate. First, we need the variance in the data set:

Sum[(Part[datasetreduced, i, 2]  -   f[
                Part[datasetreduced, i, 1]])^2, {i, 1, num}]/
      num
σ  =  Sqrt[%]

sigmaplot1  =  
      Plot[{f[t] + σ, f[t] - σ}, {t, 0, num}, 
        Frame -> True, 
        FrameLabel -> {"t", "y"}]
sigmaplot2  =  
      Plot[{f[t] + 2  σ, f[t] - 2  σ}, {t, 0, num}, 
        Frame -> True, 
        FrameLabel -> {"t", "y"}]

Show[sigmaplot1, datasetreducedplot]
Show[sigmaplot2, datasetreducedplot]

f[25]

Outlier

Here is a data set which contains an outlier:

datasetoutlier  =  {{1, 7.18437}, {2, 7.41171}, {3, 6.46758}, {4, 
            9.30980}, {5, 7.54707}, {6, 12.6665}, {7, 12.63065}, {8, 
            11.73464}, {9, 14.28273}, {10, 15.95493}, {11, 16.86929}, {12, 
            15.17158}, {13, 19.68491}, {14, 18.75987}, {15, 20.21733}, {16, 
            42.45007}, {17, 21.67026}, {18, 24.78354}, {19, 25.03960}, {20, 
            25.04889}};

datasetoutlierplot  = ListPlot[datasetoutlier]

sum  =  num
sumt  =  
      Sum[Part[datasetoutlier, i, 1], {i, 1, num}]
sumy  =  
      Sum[Part[datasetoutlier, i, 2], {i, 1, num}]
sumtt  =  
      Sum[Part[datasetoutlier, i, 1]^2, {i, 1, 
          num}]
sumty  =  
      Sum[Part[datasetoutlier, i, 1]*Part[datasetoutlier, i, 2], {i, 1, 
          num}]
delta  =  sum*sumtt  -  sumt^2

Create the linear regression model:

a  =  (sumtt*sumy - sumt*sumty)/
          delta;
b  =  (sum*sumty - sumt*sumy)/
          delta;
f[t_]  =  a  +  b*t
model  =  Plot[f[t], {t, 0, num}]

compplot  =  
      Show[model, datasetoutlierplot, Frame -> True, 
        FrameLabel -> {"t", "y"}]

ybar  =  
      Sum[Part[datasetoutlier, i, 2], {i, 1, num}]/
        num
tbar  =  
      Sum[Part[datasetoutlier, i, 1], {i, 1, num}]/
        num
Rsquare  =  
      Sum[(f[Part[datasetoutlier, i, 1]]  -   ybar)^2, {i, 1, num}]/
        Sum[(Part[datasetoutlier, i, 2] -   ybar)^2, {i, 1, 
            num}]

We see that (R^2)has been significantly reduced, to the point where we would question whether our model actually represents the data.

We can identify any outliers by computing σ and looking for data points which are more that 2σ away from the model function

Sum[(Part[datasetoutlier, i, 2]  -   f[
                  Part[datasetoutlier, i, 1]])^2, {i, 1, num}]/
        num;
σ  =  Sqrt[%]
Table[ {i, 
          Abs[Part[datasetoutlier, i, 2]  -   f[
                  Part[datasetoutlier, i, 1]] ] - 2  σ}, {i, 1, 
          num}] ;
TableForm[%]

Mathematica Package

Here is the linear regression package that comes with Mathematica. It does not encorporate the variance in each data point into the result.

<< Statistics`LinearRegression`

func[t_] = Fit[datasetreduced, {1, t}, t]

regress = 
      Regress[datasetreduced, {1, t}, t, 
        RegressionReport -> {RSquared, EstimatedVariance}]

Also, the regression package can help you identify outliers using certain options.

regress = 
      Regress[datasetoutlier, {1, t}, t, 
        RegressionReport -> {RSquared,  AdjustedRSquared, CookD, 
            PredictedResponseDelta}]

Other Linear Regressions

data  = {{0, 0.8473611502101435`}, {1, 0.16105078349847718`}, {2, 
          0.20997883486648275`}, {3, 0.5767222628194963`}, {4, 
          0.5330677559974459`}, {5, 0.7952707369812955`}, {6, 
          1.3787728263000383`}, {7, 2.0967707251295673`}, {8, 
          2.4469367135495617`}, {9, 2.7132851501919975`}, {10, 
          3.2699381172959363`}, {11, 4.595205603363667`}, {12, 
          5.268725083591067`}, {13, 5.52998002589335`}, {14, 
          6.693769672612626`}, {15, 7.735026034961229`}, {16, 
          8.039438338264981`}, {17, 9.02567289313276`}, {18, 
          10.644595538713753`}, {19, 11.48472079795253`}, {20, 
          12.540262324455885`}}

dataplot  =  
      ListPlot[data, Frame -> True, 
        FrameLabel -> {"t", "y"}]

sumtttt  =  
      Sum[Part[data, i, 1]^4, {i, 1, num}]
sumtty  =  
      Sum[Part[data, i, 1]^2*Part[data, i, 2], {i, 1, num}]

Create the linear regression model:

a  =  sumtty/sumtttt
f[t_]  =  a t^2
model  =  Plot[f[t], {t, 0, num}]

Show[dataplot, model]

Sum[(Part[data, i, 2]  -   f[Part[data, i, 1]])^2, {i, 1, 
            num}]/num;
σ  =  Sqrt[%]

sigmaplot1  =  
      Plot[{f[t] + σ, f[t] - σ}, {t, 0, num}, 
        Frame -> True, 
        FrameLabel -> {"t", "y"}]
sigmaplot2  =  
      Plot[{f[t] + 2  σ, f[t] - 2  σ}, {t, 0, num}, 
        Frame -> True, 
        FrameLabel -> {"t", "y"}]
Show[sigmaplot1, dataplot]
Show[sigmaplot2, dataplot]

regress = 
      Regress[data, {1, t^2}, t, 
        RegressionReport -> {RSquared, EstimatedVariance}]

Even More Complicated Linear Regressions

data  = {{0, 0.9043053159257124`}, {1, 8.097719084832498`}, {2, 
          6.810410294461815`}, {3, 4.632220192987015`}, {4, 
          1.1930723603688511`}, {5, -6.9574464964662495`}, {6, -5.587182715575544`}, {7, -8.665349845871182`}, {8, -5.073494692438124`}, {9, -2.3556203180876514`}, {10, 
          2.0381891304623627`}, {11, 8.480559042670984`}, {12, 
          6.550866749461399`}, {13, 
          6.065564057766881`}, {14, -0.6241447915657803`}, {15, -3.1268931737277605`}, {16, -5.2256818877843525`}, {17, -7.323617407348574`}, {18, -5.469564117752752`}, {19, -2.1255226709689263`}, {20, 
          4.496401897197252`}}
dataplot  =  
      ListPlot[data, Frame -> True, 
        FrameLabel -> {"t", "y"}]

func[t_] = 
      Fit[data, {Cos[0.2  Pi t], Sin[0.2  Pi t]}, 
        t]
funcplot  =  Plot[func[t], {t, 0, 20}]

Show[dataplot, funcplot]

regress = 
      Regress[data, {1, Cos[0.2  Pi t],  Sin[0.2  Pi t]}, t, 
        RegressionReport -> {RSquared, EstimatedVariance}]

data  =  
      Table[ {i, 3/100  i^2 + Random[Real, {0, 1}]}, {i, 0, 
          20}]
ListPlot[data]