Linear Regression

Michael Goerz
13.09.2005

Preamble

Options

SetDirectory["/home/goerz/"];

Off[General::"spell1"];

<<  Statistics`DataManipulation`

<<  Graphics`Graphics`

<< Statistics`LinearRegression`

Some Regular Expressions (As Reference)

The following shell commands can be useful for filtering datafiles:

grep '[.0-9]\|^#\|^[[:space:]]$' data.csv
This leaves any lines with data, gnuplot comment lines, and blank lines

grep '[.0-9]\|^[[:space:]]$' data.csv
This leaves any lines with data, and blank lines

grep '[.0-9]' data.csv
This leaves only lines with data

Initialization

Data is loaded from the file 'data.csv'. The file ist assumed to contain three collumns, the x-value, the y-value, and the y-error

dataset = ReadList["data.csv", {Real, Real,  Real}]

Functions

Standard Functions for Linear Regression

LinReg[daten_]:= 
Module[{x, y, Δy,  
            S , Δa0, Δb0,  a0,  b0, n, 
            Rsqared0, χSquared0},
n = Length[daten];
x = Column[daten, 1];
y = Column[daten,  2];
Δy = Column[daten, 3];
S = (∑_(i = 1)^n) 
                        (1/(((Δy[[i]]))^2)) (∑_(i = 1)^n)((((x[[i]]))^2)/(((Δy[[i]]))^2)) - ((((∑_(i = 1)^n) x[[i]]/(((Δy[[i]]))^2)))^2);
a0 = (1/S) ((∑_(i = 1)^n)((((x[[i]]))^2)/(((Δy[[i]]))^2)) (∑_(i = 1)^n) y[[i]]/(((Δy[[i]]))^2) - (∑_(i = 1)^n) 
                              x[[
                              i]]/(((Δy[[i]]))^2) (∑_(i = 1)^n)(( x[[i]] y[[i]])/(((Δy[[i]]))^2)));b0=(1/S) ((∑_(i = 1)^n) 
                            (1/(((Δy[[i]]))^2)) (∑_(i = 1)^n)(( x[[i]] y[[i]])/(((Δy[[i]]))^2)) - (∑_(i = 1)^n) 
                              x[[
                              i]]/(((Δy[[i]]))^2) (∑_(i = 1)^n)(( y[[i]])/(((Δy[[i]]))^2)));
Δa0 = (Sqrt[(1/S) (∑_(i = 1)^n)((((x[[i]]))^2)/(((Δy[[i]]))^2))]);
Δb0 = (Sqrt[(1/S) (∑_(i = 1)^n) (1/(((Δy[[i]]))^2))]);
Rsqared0 = 
              1 - ((∑_(i = 1)^n)((((y[[i]] - a0 - b0 x[[i]]))^2)/(((Δy[[i]]))^2)))/((∑_(i = 1)^n)((((y[[i]] - ((∑_(i = 1)^n) y[[i]]/(((Δy[[i]]))^2)/(∑_(i = 1)^n) (1/(((Δy[[i]]))^2)))))^2)/(((Δy[[i]]))^2)));
χSqared0 = (∑_(i = 1)^n)((((1/Δ)y[[i]] (y[[i]] - a0 - b0 x[[i]])))^2);

Print["Errors from dataset"];

Print["Number of data points: ", 
              n];

Print[
              
              TableForm[{{"", "Value",  "Error"}, {"Slope \
b",  b0,  Δb0}, {"y-Intercept a",  
                    a0,  Δa0}}]];
Print[""];
\(Print["Equation of best fit: y = ", \ b0, \ " * x + ", \ a0]\);
Print["Equation of top worst fit: y = \
",  b0 + Δb0,  " * x + ",  a0 + Δa0];
Print["Equation of bottom worst \
fit: y = ",  b0 - Δb0,  " * x + ",  
              a0 - Δa0];

\(Print["
", \*""Certainty Measure (R^2): "", \ Rsqared0]\);
\(Print["
", \*""Quality of the Errors (s^2)=((χ^2)/(n - 2)): "", χSqared0\/\(n - 2\)]\);
  

{b ->  
                b0,  Δb  ->  Δb0,  
              a ->  
                a0,  Δa ->  Δa0,  
              Rsqared ->  
                Rsqared0,  χSqared  ->  χSqared0,  
              sSquared  -> χ(Sqared0/(n - 2)) }

]

Linear Regression with Error Taken from Weighted Variance

LinRegWeightedVariance[daten_]:= 
Module[{x, y, f,  S , Δa0, Δb0,  
            a0,  b0, n, Rsqared0,  χSqared0},
n = Length[daten];
x = Column[daten, 1];
y = Column[daten,  2];
  
f = Column[daten, 3];
S = (∑_(i = 1)^n) 
                        (1/(((f[[i]]))^2)) (∑_(i = 1)^n)((((x[[i]]))^2)/(((f[[i]]))^2)) - ((((∑_(i = 1)^n) x[[i]]/(((f[[i]]))^2)))^2);
a0 = (1/S) ((∑_(i = 1)^n)((((x[[i]]))^2)/(((f[[i]]))^2)) (∑_(i = 1)^n) y[[i]]/(((f[[i]]))^2) - (∑_(i = 1)^n) x[[
                              i]]/(((f[[i]]))^2) (∑_(i = 1)^n)(( x[[i]] y[[i]])/(((f[[i]]))^2)));b0=(1/S) ((∑_(i = 1)^n) 
                            (1/(((f[[i]]))^2)) (∑_(i = 1)^n)(( x[[i]] y[[i]])/(((f[[i]]))^2)) - (∑_(i = 1)^n) x[[
                              i]]/(((f[[i]]))^2) (∑_(i = 1)^n)(( y[[i]])/(((f[[i]]))^2)));
Δa0 = (Sqrt[(1/S) (∑_(i = 1)^n)(((((x[[i]]))^2)/(((f[[i]]))^2)) (1/(n - 2)) (∑_(i = 1)^n)((((y[[i]] - a0 - b0 x[[i]]))^2)/(((f[[i]]))^2)))]);
Δb0 = (Sqrt[(1/S) (∑_(i = 1)^n)((1/(((f[[i]]))^2)) (1/(n - 2)) (∑_(i = 1)^n)((((y[[i]] - a0 - b0 x[[i]]))^2)/(((f[[i]]))^2)))]);
Rsqared0 = 
              1 - ((∑_(i = 1)^n)((((y[[i]] - a0 - b0 x[[i]]))^2)/(((f[[i]]))^2)))/((∑_(i = 1)^n)((((y[[i]] - ((∑_(i = 1)^n) y[[i]]/(((f[[i]]))^2)/(∑_(i = 1)^n) (1/(((f[[i]]))^2)))))^2)/(((f[[i]]))^2)));
χSqared0 = (∑_(i = 1)^n)((((1/f)[[i]] (y[[i]] - a0 - b0 x[[i]])))^2);
Print["Errors from variance"];

Print["Number of data points: \
", n];

Print[
              TableForm[{{"", "Value",  "Error"}, {"Slope \
b",  b0,  Δb0}, {"y-Intercept a",  
                    a0,  Δa0}}]];
Print[""];
\(Print["Equation of best fit: y = ", \ b0, \ " * x + ", \ a0]\);Print["Equation of top worst fit: y = ",  
              b0 + Δb0,  " * x + ",  
              a0 + Δa0];
Print["Equation of bottom worst fit: y \
= ",  b0 - Δb0,  " * x + ",  a0 - Δa0];
\(Print["
", \*""Quality of the Errors (s^2)=((χ^2)/(n - 2)): "", χSqared0\/\(n - 2\)]\);\(Print["
", \*""Certainty Measure \!\(R\^2\): "", \ 
              Rsqared0]\);
  

{b ->  
                b0,  Δb  ->  Δb0,  
              a ->  
                a0,  Δa ->  Δa0,  
              Rsqared ->  
                Rsqared0,  χSqared  ->  χSqared0, 
              sSquared  -> χ(Sqared0/(n - 2))}

]

Linear Regression for a Straight Line from the Origin

\(LinRegOrigin[
        daten_] := 
Module[{x, y, Δy, 
          S, Δb0, b0, n, ybar0, 
          Rsq0, χsq0}, 
n = 
          Length[daten]; 
x = 
          Column[daten, 1]; 
y = 
          Column[daten, 2]; 
Δy = 
          Column[daten, 
            3]; 
S = ∑\+\(i = 1\)\%n 1\/\((Δy[[i]])\)\^2*\(∑\+\(i = 1\)\%n\((x[[i]])\)\^2\/\((Δy[[i]])\)\^2\) - \((∑\+\(i = 1\)\%n x[[i]]\/\((Δy[[i]])\)\^2)\)\^2; 
b0 = ∑\+\(i = 1\)\%n\(\( x[[
                      i]]*
                    y[[
                      i]]\)\/\((Δy[[i]])\)\^2\)/\(∑\+\(i = 1\)\%n\((x[[i]])\)\^2\/\((Δy[[i]])\)\^2\)\ ; 
Δb0 = 
          Sqrt[1/\(∑\+\(i = 1\)\%n\((x[[i]])\)\^2\/\((Δy[[i]])\)\^2\)]; 
ybar0 = ∑\+\(i = 1\)\%n\( y[[
                    i]]\/\((Δy[[i]])\)\^2\)/\(∑\+\(i = 1\)\%n 1\/\((Δy[[i]])\)\^2\); 
Rsq0 = 
          1 - ∑\+\(i = 1\)\%n\(\((y[[i]] - b0*x[[i]])\)\^2\/\((Δy[[i]])\)\^2\)/\(∑\+\(i = 1\)\%n\((y[[i]] - ybar0)\)\^2\/\((Δy[[i]])\)\^2\); 
χsq0 = ∑\+\(i = 1\)\%n\((y[[i]] - b0*x[[i]])\)\^2\/\((Δy[[i]])\)\^2; 
Print["Errors from dataset"]; 
Print["Number of data points: ", 
          n]; 
Print[
          TableForm[{{"", "Value", "Error"}, {"Slope b", 
                b0, Δb0}}]]; 
Print[""]; 
Print["Equation of best fit: y = ", \ 
          b0, \ " * x"]; 
Print["Equation of top worst fit: y = ", \ 
          b0 + Δb0, \ " * x"]; 
Print["Equation of bottom worst fit: y = ", \ 
          b0 - Δb0, \ " * x"]; 
Print[""]; 
Print[\*""Certainty Measure (R^2): "", 
          Rsq0]; 
Print[\*""Quality of Errors (s^2)=((χ^2)/(n - 1)): "", χsq0\/\(n - 1\)]; \ 
{a -> 0, Δa -> \ 0, 
          b -> \ b0, \ Δb -> Δb0, 
          Rsqared -> Rsq0, \ χSqared\  -> \ χsq0, 
          sSquared\  -> χsq0\/\(n - 2\)}
]\)

Instructions for Exponential Data

It is best to linearize exponential data first with the following command:

linearize[expdata_] := 
        Replace[expdata, {xi_, yi_, Δyi_} ->  {xi, 
              Log[yi], (1/yi)*Δyi}, 2];

Then, you can apply the functions for linear data. If you want to plot the data as exponential functions later, simply raise the full right hand side of the fit equation into the exponent.

Exporting Data

To export data to a file, use
Export["data.csv", dataset,  "Table"]

Application

Calculation

parameter = LinReg[dataset];

Plot

\(\(DisplayTogether[
ErrorListPlot[
          dataset], 
Plot[{b\ x + a} /. parameter // 
            Evaluate, {x, 0, Max[Column[dataset, 1]]}, 
          PlotStyle -> {Dashing[{0}]}], 
        Plot[{\((b + Δb)\) 
                    x + \((a + Δa)\), \((b - Δb)\) x + \((a - Δa)\)} /. parameter // 
            Evaluate, {x, 0, 
            Max[Column[dataset, 
                1]]}, 
PlotStyle -> {Dashing[{0.01}]}], ImageSize -> 700, PlotLabel -> "dataset", 
        AxesLabel -> {"x", "y"}];\)\)