## [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. }