Python: Least Square Adjustment; Example 5

20110419

find parameter of function y; f(y)

f(y) = a[0] + a[1]*x + a[2]*x*x + a[3]*x*x*x

data is
x = 1, 2, 3, 4, 5  and x is constance
y = 1.74, 2.79, 4.33, 6.16, 8.51

n=5, n[0]=4, r=1

condition equation

y[i] + v[i] = a[0] + a[1]*x[i] + a[2]*x[i]*x[i] + a[3]*x[i]*x[i]*x[i]

then

v[i] = a[0] + a[1]*x[i] + a[2]*x[i]*x[i] + a[3]*x[i]*x[i]*x[i] – y[i]

Solve by Python

## Least Square Adjustment
## Example: 5
## Adjustment by Parameter
## Leveling

''' data x,y fix x
x = 1., 2., 3., 4., 5.
y = 1.74, 2.79, 4.33, 6.16, 8.51

f(y) = a_0 + a_1*x + a_2*x^2 + a_3*x^3
'''

## n= 5, n_0=4, r = n-n_0 = 1
''' condition equation
v_1 = a_0 + a_1*x_1 + a_2*x_1^2 + a_3*x_1^3 - y_1
v_2 = a_0 + a_1*x_2 + a_2*x_2^2 + a_3*x_2^3 - y_2
v_3 = a_0 + a_1*x_3 + a_2*x_3^2 + a_3*x_3^3 - y_3
v_4 = a_0 + a_1*x_4 + a_2*x_4^2 + a_3*x_4^3 - y_4
v_5 = a_0 + a_1*x_5 + a_2*x_5^2 + a_3*x_5^3 - y_5'''

import numpy as np

##set data
x = (1.,2.,3.,4.,5.)
y = (1.74, 2.79, 4.33, 6.16, 8.51)

##general equation for v
##v[i] = a_0 + a_1*x[i] + a_2*x[i]^2 + a_3*x[i]^3 - y[i]
## set data to v + B * Delta = f
## set all a = 1 for find index of a
a = (1.,1.,1.,1.)

B = []
for i in range(5):
 B.append([-a[0],-a[1]*x[i],-a[2]*x[i]*x[i],-a[3]*x[i]*x[i]*x[i]])

''' from v + B * Delta = f
y[i] + v[i] - a[0] - a[1]*x[i] - a[2]*x[i]*x[i] - a[3]*x[i]*x[i]*x[i] = 0
v[i] - a[0] - a[1]*x[i] - a[2]*x[i]*x[i] - a[3]*x[i]*x[i]*x[i] = -y[i]
v[i] + [-a[0] - a[1]*x[i] - a[2]*x[i]*x[i] - a[3]*x[i]*x[i]*x[i]] = -y[i]
'''
f = np.asmatrix(np.dot(y,-1)).reshape(len(y),1)

## weight matrix
W = np.matrix([[1./np.square(0.02),0.,0.,0.,0.],[0,1./np.square(0.02),0,0,0],[0,0,1./np.square(0.02),0,0],[0,0,0,1./np.square(0.04),0],[0,0,0,0,1./np.square(0.04)]])

B = np.asmatrix(B)
## N = B^T * W * B
N = B.T * W * B

## t = B^T * W * f
t = B.T * W * f

## delta = N^(-1) * t
invN = np.linalg.inv(N)
delta = invN * t

print "parameter"
print delta

result is
parameter
[[ 1.13834711]
 [ 0.35227273]
 [ 0.25132231]
 [-0.00549587]]

parameter is a[0], a[1], a[2] and a[3]

parameter
a[0] = 1.13834711
a[1] = 0.35227273
a[2] = 0.25132231
a[3] = -0.00549587
Ans

Advertisement