[REBOL] Re: GPS: minimizing sum of squares algorithm
From: jan:skibinski:sympatico:ca at: 10-Nov-2002 3:40
Jose,
As promised, here is a solution to your GPS minimization
problem, or rather a part of it. The only function I am
attaching here is the main GPS function. I have not tested
it, so it might be buggy, therefore threat it as the outline
only. I developed and tested today various pieces that
the gps function is calling, including a solver for a system
of linear equations. If you want to use them just ask for
it and I will send it to you. I did not plan to develop any
linear algebra package at this time and whatever I developed
today is not optimized and quite rough around the edges.
For this reason I am not publicly posting it to the library.
But you can use any other solver, or any other package with
the gps function below.
All the best,
Jan
comment {
Georeference GPS data modeling:
------------------------------
Find six constants relating [x y] coordinates
and [long lat] coordinates, which we abbreviate
here as [u v]:
x = a1 + a2 u + a3 v
y = b1 + b2 u + b3 v
To find a1, a2, a3, b1, b2, b3 we need at least 3 sets
of data for x, y, u and v - the more the merrier:
[x1 x2 x3 ... xn]
[y1 y2 y3 ... yn]
[u1 u2 u3 ... un]
[v1 v2 v3 ... vn]
Minimization of the least square error F:
F = Sum (xi - a1 - a2 ui - a3 vi)**2
+ (yi - b1 - b2 ui - b3 vi) ** 2
Six stationary conditions for extremum of F:
dF/da1 = 0, dF/da2=0, dF/a3=0,
dF/db1=0, dF/db2=0, dF/db3=0
Simplification:
-------------
I assumed below the same average variance for all the
samples,which is good enough if your data is reasonably
correct, otherwise the model would need to deal with
individual variances, and the description would not be
as simple as this one.
And because the variance is constant it will not appear
anywhere in the equations that follow.
Since variables x and y are independent the system of
six equation separates itself into two system of three
equation each.
We have to separately solve two matrix equations,
M a = p
M b = q
The matrix M is:
| n | Sum ui | Sum vi |
| | | |
| Sum ui | Sum ui ui | Sum ui vi |
| | | |
| Sum vi | Sum ui vi | Sum vi vi |
Vector p-transposed = [Sum xi, Sum xi ui, Sum xi vi]
Vector q-transposed = [Sum yi, Sum yi ui, Sum yi vi]
Internally, we will represent matrix M as block of
three block-columns.
}
gps: func [
{A block of two blocks [[a1 a2 a3][b1 b2 b3]]
containing estimated model parameters, computed
from 4 sets of input data 'xs 'ys 'us 'vs
}
xs [block!] {[x1 x2 .... xn]}
ys [block!] {[y1 y2 ..... yn]}
us [block!] {[u1 u2 ..... un]}
vs [block!] {[v1 v2 ..... vn]}
/local n m p q result [block!]
][
simple-sum: func[x][foldl :+ 0 x]
n: length? xs
m: copy []
append/only m reduce [n simple-sum us simple-sum vs]
append/only m reduce [simple-sum us product-vv us us product-vv
us vs]
append/only m reduce [simple-sum vs product-vv us vs product-vv
vs vs]
p: reduce [simple-sum xs product-vv xs us product-vv xs vs]
q: reduce [simple-sum ys product-vv ys us product-vv ys vs]
result: copy []
append/only result solve m p
append/only result solve m q
result
]
comment {
'product-vv stands for scalar inner product of two vectors.
'Solve uses my home grown algorithm based on Gram-Schmidt
orthogonalization procedure. A typical approach here would be
to use the LU decomposition technique with pivoting.
}