Curve Fitting

Example:

Given a set of points (x1, y1), …….(xn, yn) we want to fit these points to an exponential curve, such that

F(x) = C eAx

Using least square procedure

E(C, A) = S ( f(xk) – Yk ) 2 where k = 1 to n

Since f(x) = C eAx , we have

E(C, A) = S ( C eAxk – Yk ) 2

Then we take the derivative with respect to C

cE/cC = å 2 (CeAxk -y) eAxk

= 2 å (Ce2Axk - eAxk yk) = 0

= C å (e2Axk) - å (eAxkyk) = 0…………….1

Then we take another derivative with respect to A

cE/cA =å 2 (CeAxk -y) CxkeAxk = 0

= å (C2xke2Axk-ykCxkeAxk) = 0

= C2å(xke2Axk) – C å(ykxkeAxk) = 0………2

1 and 2 give us a system of two non-linear system equation can ve solved using Newton method.

Using Linearization

Basic idea: transform the equation from y = C eAx into y = A ln(x) + B

Exp:

y = C eAx ……………… put both sides into ln form, then we have

Ln(y) = ln( C eAx )

Ln(y) = ln(C) + Ax

Then we have Z = ln(y) , ln(c) = B, transform it into

Z = Ax + B

Exp

Given data, fit these into y = C eAx

xk         yk

1.5

2.5

3.5

5.0

7.5

 

Solution:

Xk Yk Zk = ln(y) (Xk)2 XkZk

0 1.5 .4055 0 0

1 2.5 .9163 1 .9163

2 3.5 1.2528 4 2.5056

3 5.0 1.6094 9 4.8282

4 7.5 2.0149 16 4.0599

S Xk = 10

S Zk = 6.1989

S (Xk)2 = 30

S XkZk = 16.3097

Plugging in all those values into 1 and 2, we have

0 = A ( 30 ) + B ( 10 ) – 16.3097 ………….1

0 = A ( 10 ) + B ( 5 ) – 6.1989 …………….2

Solve those equations simultaneously, then we have A = .3913 and B = .4574

Since C = eB , C = 1.579

Our equation is y = .3913 e 1.579x

Linear Least Square

A set of N data points (x1,y1), …, (xn, yn), we want to:

f(x) = c1f1(x) + c2f2(x) + …+ cMfM(x)

We want c1, c2, …, cM to minimize the square error to approximate (x1, y1), …,(xn,yn)

 E(c1, c2, …, cM) = S( f(xk)-yk )2 where k is 1 to n

= S(c1f1(x k) + c2f2(x k) + …+ cMfM(x k) – y k) 2 where k = 1 to n

= S S (c jf j(x k) – y k) 2 where k = 1 to n and j = 1 to m

Apply the same routine we take the derivative with respect to ci where i = 1 to n and j = 1 to m

E/ci = S 2*S(c jf j(x k) – y k)*(fi (x k)) = 0

= S S(c jf j(xk) fi (x k)) – S (y kfi (x k)) = 0

Solve this systems we get c1, c2, …, cM to fit (x1, y1), …,(xN,yN)

f(x) = c1f1(x) + c2f2(x) + …+ cMfM(x)

Polynomial fitting

Given a set of data points (x1, y1), …,(xN,yN), we want to find the c1, c2, …, cM+1 such that

f(x) = c1+ c2x + c3x2 + …+ cM+1xM

Example:

We want to find A, B, C that satisfy

f(x) = Ax2 + Bx + C

that best fits (x1, y1), …,(xn,yn) minimizing the square error. 

E(A,B,C) = S (Ax k 2 + Bx k + C – y k) 2 where k = 1 to n

 

Take the derivative with respect of A

E/A = S 2*(Ax k 2 + Bx k + C – y k)* xk 2 = 0

AS xk4 + BS xk3 + CS xk2 = S xk2 yk………………………..1

Take the derivative with respect of B

E/B = S 2*(Ax k 2 + Bx k + C – y k)* xk = 0

ASk=1N xk3 + BSk=1N xk2 + CSk=1N xk = Sk=1N xk yk…………2

Take the derivative with respect of A

E/C = S 2*(Ax k 2 + Bx k + C – y k) = 0

ASxk2 + BS xk + CN = S yk………………………………..3

1, 2, and 3 form a system of 3 linear equations of 3 variables which enable us to find our solution

f(x) = Ax2 + Bx + C

The problem with using high degree polynomials to approximate a set of data is that the higher the degree of the polynomial, the more maximum and minimum points it will have, making the curve wiggle.  To avoid this problem, we use spline functions. 

Spline functions

Spline are a set of curves that fit a set of points piecewise. Instead of having one single polynomial to fit all the points we have multiple polynomials each with certain intervals.

Piecewise Cubic Splines

To have continuity also in the first and second derivatives, we use points before and beyond the interval of the spline.

Properties of cubic splines:

S(Xk) = yk, guarantee that the splines pass through the data points

Sk(Xk+1) = Sk+1(Xk+1) , continuous in y

Sk’ (Xk+1) = Sk+1’(Xk+1) , continuity in first derivatives

Sk’’(Xk+1) = Sk+1’’(Xk+1) , continuity in second derivatives

Suppose we want to fit N+1 data points (x0, y0), …, (xN, yN), where a = x0 < x1 < … < xn = b, into N cubic polynomials Sk(x) with coefficients Sk,0, Sk,1, Sk,2 and Sk,3

We want Sk(x) k = 0, …, N-1 to satisfy

I.                   S(x) = Sk(x) = Sk,0 + Sk,1(x-xk) + Sk,2(x- xk)2 + Sk,3 (x- xk)3

for xk £ x £ xk+1 and k = 0, 1, …, N-1

II.                 S(xk) = yk  for k = 0, 1, …, N

III.               Sk(xk+1) = Sk+1(xk+1)  for k = 0, 1, …, N-2

IV.              S’k(xk+1) = S’k+1(xk+1)  for k = 0, 1, …, N-2

V.                 S"k(xk+1) = S"k+1(xk+1)  for k = 0, 1, …, N-2 

For each polynomial Sk(x) there are four coefficients to be determined. 

Sk,0, Sk,1, Sk,2, Sk,3 (4 degrees of freedom)

Also we have N splines S0(x), S1(x), …,  SN(x)

Since Sk(x) is cubic, the second derivative is a line.  We can obtain the equation of this line using Lagrange approximation. 

S"k(x) = S"(xk)(x-xk+1)/(xk- xk+1) + S"(xk+1)(x- xk)/(xk+1- xk)

So we have,

S"k(x) = S"(xk) when x = xk

S"k(x) = S"(xk+1) when x = xk+1

Assume: (we create constants)

mk = S"(xk)

mk+1 = S"(xk+1)

hk = xk+1- xk

S"k(x) = mk(x- xk+1)/- hk  + mk+1(x- xk)/ hk

= mk(xk+1-x)/hk  + mk+1(x- xk)/ hk

To obtain Sk(x) we integrate twice:

Sk(x) = òò S"(x) dx dx = òò mk(xk+1-x)/hk  + mk+1(x- xk)/ hk dx dx

= ò -½(mk / hk)( xk+1-x)2 + ½(mk+1/hk)(x-xk)2 + C1 dx

= ½(1/3)(mk/hk)(xk+1-x)3 + ½(1/3)(mk+1/hk)(x- xk)3 + C1x + C2

= (mk/6hk)(xk+1-x)3 + (mk+1/6hk)(x-xk)3 + pk(xk+1-x) + qk(x- xk)

Substitute xk, xk+1 and S(xk) = yk and S(xk+1) = yk+1

Sk(xk) = yk = (mk/6hk)(xk+1-xk)3 + (mk+1/6hk)(xk-xk)3 + pk(xk+1-xk) + qk(xk- xk)

yk = (mk/6hk)(hk)3 + pk(hk) = 1/6(mkhk2) + pkhk 

Sk(xk+1) = yk+1 = (mk/6hk)(xk+1-xk+1)3 + (mk+1/6hk)(xk+1-xk)3 + pk(xk+1-xk+1) + qk(xk+1- xk)

= (mk+1/6hk)(hk)3 + qk(hk) = 1/6(mk+1hk2) + qkhk

Yielding us to

yk = 1/6(mkhk2) + pkhk………………………………..1

yk+1 = 1/6(mk+1hk2) + qkhk……………………………2

We use 1 to obtain pk

pk = (yk- mkhk2/6)(1/ hk) = (yk/ hk- mkhk/6)

We use 2 to obtain qk

qk = (yk+1-mk+1hk2/6)(1/hk) = (yk+1/hk -mk+1hk/6)

Substitute 1 and 2 into

Sk(x) = mk/6hk)(xk+1-x)3 + (mk+1/6hk)(x-xk)3 + pk(xk+1-x) + qk(x- xk)

Yields

Sk(x) = (mk/6hk)(xk+1-x)3 + (mk+1/6hk)(x-xk)3 + (yk/hk- mkhk/6)(xk+1-x) + (yk+1/hk-mk+1hk/6)(x- xk)