Understanding Least Squares and Data Modeling

Data Modeling Patrice Koehl Department of Biological Sciences National University of Singapore http://www.cs.ucdavis.edu/~koehl/Teaching/BL5229 koehl@cs.ucdavis.edu

Data Modeling Data Modeling: least squares Data Modeling: Non linear least squares Data Modeling: robust estimation

Data Modeling Data Modeling: least squares Linear least squares

Least squares Suppose that we are fitting N data points (xi,yi) (with errors i on each data point) to a model Y defined with M parameters aj: Y(x;a1,a2,...,aM) The standard procedure is least squares: the fitted values for the parameters aj are those that minimize: 2 N yi-Y(x;a1,...,aM) si c2= i=1 Where does this come from?

Least squares Let us suppose that: The data points are independent of each other Each data point has a measurement error that is random, distributed as a Gaussian distribution around the true value Y(xi) The probability of the data points, given the model Y is then: 2 N yi-Y(xi) si exp -1 P(data/Model) 2 i=1

Least squares Application of Bayes s theorem: P(Model/Data) P(Data/Model)P(Model) With no information on the models, we can assume that the prior probability P(Model) is constant. Finding the coefficients a1, aM that maximizes P(Model/Data) is then equivalent to finding the coefficients that maximizes P(Data/Model). This is equivalent to maximizing its logarithm, or minimizing the negative of its logarithm, namely: 2 N yi-Y(x) si 1 2 i=1

Fitting data to a straight line

Fitting data to a straight line This is the simplest case: Y(x) =ax+b Then: 2 N yi-axi-b si c2= i=1 The parameters a and b are obtained from the two equations: xiyi-axi-b ( si yi-axi-b si ) N c2 a c2 b =0 = -2 2 i=1 N =0 = -2 2 i=1

Fitting data to a straight line Let us define: N N N N N 2 1 xi si yi si xi si xiyi si S = Sx= Sy= Sxx= Sxy= si 2 2 2 2 2 i=1 i=1 i=1 i=1 i=1 aSxx+bSx aSx+bS = = Sxy Sy then a =SxyS -SxSy SxxS -Sx b =SxxSy-SxSxy SxxS -Sx a and b are given by: 2 2

Fitting data to a straight line We are not done! S 2= sa SSxx- Sx Sx SSxx- Sx 2 Uncertainty on the values of a and b: sb 2= 2 Evaluate goodness of fit: -Compute 2 and compare to N-M (here N-2) -Compute residual error on each data point: Y(xi)-yi -Compute correlation coefficient R2

Fitting data to a straight line

General Least Squares Y(x) =a1X1(x)+a2X2(x)+...+aMXM(x) Then: 2 N yi-a1X1(xi)-...-aMXM(xi) si c2= i=1 The minimization of 2 occurs when the derivatives of 2 with respect to the parameters a1, aMare 0. This leads to M equations: N c2 ak 1 ( )Xk(xi) =0 = yi-a1X1(xi)-...-aMXM(xi) si i=1

General Least Squares Aij=Xj(xi) Define design matrix A such that si

General Least Squares bi=yi Define two vectors b and a such that si and a contains the parameters Note that 2 can be rewritten as: c2= Aa -b 2 The parameters a that minimize 2 satisfy: ( )a = ATb ATA These are the normal equations for the linear least square problem.

General Least Squares How to solve a general least square problem: 1) Build the design matrix A and the vector b 2) Find parameters a1, aMthat minimize c2= Aa -b 2 (usually solve the normal equations) 3) Compute uncertainty on each parameter aj: s(aj)2=C-1(j,j) if C = ATA, then

Data Modeling

Non linear least squares In the general case, g(X1, ,Xn) is a non linear function of the parameters X1, Xn; 2 is then also a non linear function of these parameters: ( ) 2 Yi-g X1,...,Xn;t si N ( )= c2= f X1,...,Xn i=1 Finding the parameters X1, ,Xn is then treated as finding X1, ,Xn that minimize 2.

Minimizing 2 Some definitions: Gradient: The gradient of a smooth function f with continuous first and second derivatives is defined as: x f f f ( ) X = ... ... f x x 1 i N Hessian The n x n symmetric matrix of second derivatives, H(x), is called the Hessian: 1 ... ... f x H x x 2 2 2 f f f ... ... ... 2 x x x 1 ... 1 j N ... x x x 2 2 2 f f = ... ... ( ) ... x x x 1 i i j i N ... ... ... ... x x 2 2 2 f f f ... ... 2 N x x x 1 N N j

Minimizing 2 Minimization of a multi-variable function is usually an iterative process, in which updates of the state variable x are computed using the gradient and in some (favorable) cases the Hessian. Steepest descent (SD): The simplest iteration scheme consists of following the steepest descent direction: ( sets the minimum along the line defined by the gradient) = - a ( ) x x f x + 1 k k k Usually, SD methods leads to improvement quickly, but then exhibit slow progress toward a solution. They are commonly recommended for initial minimization iterations, when the starting function and gradient-norm values are very large.

Minimizing 2 Conjugate gradients (CG): In each step of conjugate gradient methods, a search vector pk is defined by a recursive formula: ( ) k x f p 1 + - = + b p + 1 k k k The corresponding new position is found by line minimization along pk: x x + = +1 l p k k k k the CG methods differ in their definition of .

Minimizing 2 Newton s methods: Newton s method is a popular iterative method for finding the 0 of a one-dimensional function: ( ) ( ) x ' g x = - k x x + 1 k k g k x3 x2 x1 x0 It can be adapted to the minimization of a one dimensional function, in which case the iteration formula is: ( ) ( ) x ' ' ' g x = - k x x + 1 k k g k

Minimizing 2 The equivalent iterative scheme for multivariate functions is based on: ( ) ( ) k x - = - 1 x x H f x + 1 k k k k Several implementations of Newton s method exist, that avoid Computing the full Hessian matrix: quasi-Newton, truncated Newton, adopted-basis Newton-Raphson (ABNR),

Data analysis and Data Modeling Data Modeling: robust estimation

Robust estimation of parameters Least squares modeling assume a Gaussian statistics for the experimental data points; this may not always be true however. There are other possible distributions that may lead to better models in some cases. One of the most popular alternatives is to use a distribution of the form: r(x) =e- x Let us look again at the simple case of fitting a straight line in a set of data points (ti,Yi), which is now written as finding a and b that minimize: N Z(a,b) = Yi-ati-b i=1 b = median(Y-at) and a is found by non linear minimization

Robust estimation of parameters