Least squares fit for Linear, Square polynomial and Circular (with deduction)

The basic approach of the Least squares fit is to minimize the function f(x) by varying the parameter p:

min ∑ (y_i − f_p(x_i))²

Variation of p:

ddp ∑ (y_i − f_p(x_i))² = 0

Linear case

f(x) = kx + b

Minimization core L:

L = 1/N × ∑ [k²(x²) + 2kbx − 2k(xy) − 2by + b² + (y²)]

Arrange and get derivative:

L(b) = 2kbx̄ − 2bȳ + b²

L'(b) = 2kx̄ − 2ȳ + 2b = 0

b = ȳ − kx̄

L(k) = k²∑(x²)/N + 2kb∑x/N − 2k∑(xy)/N

L'(k) = 2k∑(x²)/N + 2bx̄ − 2∑(xy)/N = 0

S_xy = 1/N × ∑ xy

S_xx = 1/N × ∑ x²

kS_xx + bx̄ − S_xy = kS_xx + (ȳ − kx̄)x̄ − S_xy = kS_xx + ȳx̄ − kx̄² − S_xy = 0

kS_xx − kx̄² = S_xy − ȳx̄

k = S_xy − x̄ȳ S_xx − x̄²

Square polynomial (Quadratic polynomial) case

f(x) = ax² + bx + c

Minimization core Q:

Q = 1/N × ∑ [a²(x⁴) + 2ab(x³) − 2a(x²y) + b²(x²) + 2ac(x²) − 2b(xy) + 2bcx − 2cy + c² + y²]

Arrange and get derivative:

Q(c) = c² + 2ac∑x²/N + 2bc∑x/N − 2c∑y/N

S_xx = 1/N × ∑ x²

Q'(c) = 2c + 2aS_xx + 2bx̄ − 2ȳ = 0

c = ȳ − aS_xx − bx̄

Q(b) = b²∑x²/N + 2ab∑x³/N − 2b∑xy/N + 2bc∑x/N = b²S_xx + 2ab∑x³/N − 2b∑xy/N + 2b(ȳ − aS_xx − bx̄)x̄ = b²S_xx + 2ab∑x³/N − 2b∑xy/N + 2bx̄ȳ − 2abx̄S_xx − 2b²x̄²;

S_xxx = 1/N × ∑ x³

S_xy = 1/N × ∑ xy

Q'(b) = 2bS_xx + 2aS_xxx − 2S_xy + 2x̄ȳ − 2ax̄S_xx − 4bx̄² = 0

Q'(b) = bS_xx − bx̄² + aS_xxx − S_xy + x̄ȳ − ax̄S_xx = 0

b(S_xx − x̄²) = − aS_xxx + S_xy − x̄ȳ + ax̄S_xx

R₀ = (S_xy − x̄ȳ)/(S_xx − x̄²)

R₁ = (S_xxx − x̄S_xx)/(S_xx − x̄²)

b = R₀ − aR₁

Q(a) = a²x⁴ + 2abx³ - 2ax²y + 2acx² = a²x⁴ + 2abx³ - 2ax²y + 2a(ȳ - aS_xx - bx̄)x² = a²x⁴ + 2abx³ - 2ax²y + 2aȳx² - 2a²S_xxx² - 2abx̄x²

Q(a) = a²x⁴ - 2ax²y + 2aȳx² - 2a²S_xxx² + 2aR₀(x³ - x̄x²)-2a²R₁(x³ - x̄x²)

Q'(a) = 2a∑x⁴/N - 2∑x²y/N + 2ȳ∑x²/N - 4aS_xx∑x²/N + 2R₀∑(x³ - x̄x²)/N - 4aR₁∑(x³ - x̄x²)/N = 0

S_xxxx = 1/N × ∑ x⁴

S_xxy = 1/N × ∑ x²y

S₂ = S_xxx - x̄S_xx

2aS_xxxx - 2S_xxy + 2ȳS_xx - 4aS_xx² + 2R₀S₂ - 4aR₁S₂ = 0

aS_xxxx - 2aR₁S₂ - S_xxy + ȳS_xx - 2aS_xx² + R₀S₂ = 0

a(S_xxxx - 2R₁S₂ - 2S_xx²) - S_xxy + ȳS_xx + R₀S₂ = 0

a = S_xxy - ȳS_xx - R₀S₂ S_xxxx - 2R₁S₂ - 2S_xx²

Circle case

f(x) = (x−x₀)² + (y−y₀)² − R²

Minimization core C all terms (doubled above diagonal, x02 means x₀²):

	x2	-2x0x	x02	y2	-2y0y	y02	-R2
x2	x4	-2x0x3	2x02x2	2x2y2	-4y0x2y	2y02x2	-2R2x2
-2x0x		4x02x2	-4x03x	-4x0xy2	8x0y0xy	-4x0y02x	4R2x0x
x02			x04	2x02y2	-4x02y0y	2x02y02	-2R2x02
y2				y4	-4y0y3	2y02y2	-2R2y2
-2y0y					4y02y2	-4y03y	4R2y0y
y02						y04	-2R2y02
-R2							R4

C(R) = -2R²x² +4R²x₀x -2R²x₀² -2R²y² +4R²y₀y -2R²y₀² +R⁴

C(x0) = x₀⁴ - 4x₀³x + 6x₀²x² - 2x₀x³ - 4x₀xy² + 2x₀²y² + 2x₀²y₀² - 4x₀²y₀y - 4x₀y₀²x + 8x₀y₀xy - 2R²x₀² + 4R²x₀x

C(y0) = y₀⁴ - 4y³y + 6y₀²y² - 2y₀y³ - 4y₀x²y + 2y₀²y² + 2x₀²y₀² - 4x₀²y₀y - 4x₀y₀²x + 8x₀y₀xy - 2R²y₀² + 4R²y₀y

Arrange and get derivative by R

C'(R) = 1/N ∑ [-4Rx² +8Rx₀x -4Rx₀² -4Ry² +8Ry₀y -4Ry₀² +4R³] = 0

As R > 0 by define, divide by 4R:

C'(R) =1/N ∑[-x² +2x₀x -x₀² -y² +2y₀y -y₀² +R²] = 0

R² = ∑x² -2x₀∑x +x₀² +∑y² -2y₀∑y +y₀² = S_xx - 2x₀x̄ + x₀² + S_yy - 2y₀ȳ + y₀²

Mean optimization

Instead of using large sums of squares (and other powers), sum of squares about the mean can be used to ease the computational task and reduce rigidity. However, the result is obtained in two passes.

Example transfer formula for S_xx:

S_xx = 1/N × ∑ x²

x̄ = 1/N × ∑ x

SS_xx = 1/N × ∑ (x − x̄)² = 1/N × ∑ x² − 1/N × 2 × x̄ ∑ x + 1/N × x̄² ∑ 1 = 1/N × ∑ x² − 2x̄² + x̄²

SS_xx = S_xx − x̄²

S_xx = SS_xx + x̄²

Linear case

SS_xx = 1/N × ∑ (x − x̄)²

SS_xy = 1/N × ∑ (x − x̄)(y − ȳ)

SS_xy = S_xy − x̄ȳ

k = S_xy − x̄ȳ S_xx − x̄² = SS_xySS_xx

b = ȳ − k x̄

Square case

SS_xxx = 1/N × ∑ (x − x̄)³

SS_xxxx = 1/N × ∑ (x − x̄)⁴

SS_xxy = 1/N × ∑ (x − x̄)²(y − ȳ)

SS_xxx = S_xxx − 3x̄S_xx + 2x̄³ = S_xxx − x̄S_xx - 2x̄SS_xx

SS_xxxx = S_xxxx − 4x̄S_xxx + 6x̄²S_xx − 3x̄⁴

SS_xxy = S_xxy − ȳS_xx - 2x̄S_xy + 2x̄²ȳ

r₁ = SS_xxx SS_xx

r₂ = SS_xy SS_xx

Variables:

a = SS_xxy - r₁SS_xy SS_xxxx - r₁SS_xxx - (SS_xx)²

Circle case

SS_yy = 1/N × ∑ (y − ȳ)²

SS_yyy = 1/N × ∑ (y − ȳ)³

SS_xyy = 1/N × ∑ (x − x̄)(y − ȳ)²

D = SS_xxSS_yy - (SS_xy)²

s₁ = SS_xxx + 2x̄SS_xx + SS_xyy + 2ȳSS_xy

s₂ = SS_yyy + 2ȳSS_yy + SS_xxy + 2x̄SS_xy

x₀ = SS_yy s₁ - SS_xy s₂ 2D

y₀ = SS_xx s₂ - SS_xy s₁ 2D

R² = SS_xx+(x̄-x₀)² + SS_yy+(ȳ-y₀)²

R = sqrt(R²)

Errors estimation

Summation

Summation of a large number of values over a wide range can lead to the accumulation of a significant error. Use Kahan summation algorithm wiki

Tcl Code example

# Sums $ss??? are without 1/N factor
# Linear fit  y=a1*x+a0
	set a1 [expr 1.0*$ssxy/$ssxx]
	set a0 [expr 1.0*($my-$a1*$mx)]
# Simple square fit  y=a2*x^2+a0
	set a2 [expr $ssxy/($ssxxx+2.0*$mx*$ssxx)]
	set a0 [expr $my-$a2*($ssxx/$n+$mx*$mx)]
# Full square fit y=a2*x^2+a1*x+a0
	set r1 [expr $ssxxx/$ssxx]
	set r2 [expr $ssxy/$ssxx]
	set r3 [expr $ssxx/$n]
	set a2 [expr ($r1*$ssxy-$ssxxy)/($ssxx*$r3+$ssxxx*$r1-$ssxxxx)]
	set a1 [expr $r2-$a2*($r1+2.0*$mx)]
	set a0 [expr $my-$a1*$mx-$a2*($r3+$mx*$mx)]
# Circle fit (sums 1/N factored)
	set det [expr $ssxx*$ssyy-$ssxy*$ssxy]
	set s1  [expr $ssxxx+2.0*$mx*$ssxx+$ssxyy+2.0*$my*$ssxy]
	set s2  [expr $ssyyy+2.0*$my*$ssyy+$ssxxy+2.0*$mx*$ssxy]
	set x0  [expr ($ssyy*$s1-$ssxy*$s2)/(2.0*$det)]
	set y0  [expr ($ssxx*$s2-$ssxy*$s1)/(2.0*$det)]
	set rr  [expr $ssxx+$mx*$mx-2.0*$mx*$x0+$x0*$x0+$ssyy+$my*$my-2.0*$my*$y0+$y0*$y0]
	set r0  [expr sqrt($rr)]