GPS: minimizing sum of squares algorithm

Has anybody written this algorithm in Rebol ? I need this to georeference gps data, i.e. find a, b, c, d, e, f below given three calibration points x = a + b lat + c lon y = d + e lat + f lon It's available as a math library in most languages but maybe some rebol guru has it too ! Regards

Jose, I could certainly help you with this but I would need to know more about the shape of the functions b c e f in order to model them appriopriately. For example, depending on what those functions represent a polynomial approximation might be of no use or it might be quite good actually. A reference to existing georeference algorithms could be a good starting point. Regards, Jan jose wrote:

Hi Jose I implemented some linear alrerbra functions a few years back. I will see if I can dig them up over the weekend On Fri, 8 Nov 2002, [iso-8859-1] jose wrote:

a,b,c,d,e,f are constants. I basically want to solve a set of 6 equations with six unknowns, using one of the existing linear algebra algorithms. The algorithm implementation I'm looking for is also known as RMS, root mean squares but couldn't find a good algorithm description in the web Regards Jose --- Jan Skibinski <[jan--skibinski--sympatico--ca]> escribió: > Jose,

Jose,

> a,b,c,d,e,f are constants. > I basically want to solve a set of 6 equations with > Regards > Jose

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

Jan, Many thanks for your help. Can you pls send me your functions so that I can test your gps functions I find that there are very few math related scripts available in rebol, so your scripts are most valuable! Thanks --- Jan Skibinski <[jan--skibinski--sympatico--ca]> escribió: > Jose,

Jose and all,

> Many thanks for your help. Can you pls send me your > functions so that I can test your gps functions > > I find that there are very few math related scripts > available in rebol, so your scripts are most valuable! >

Can you pls send a link to your web article Indexless Linear Algebra Algorithms ? The link at numeric-quest seems broken ! Thanks again --- Jan Skibinski <[jan--skibinski--sympatico--ca]> escribió: >

Hi Jose, Hmmm I recall "root mean square" from Calculus not Linear and it was the sqrt of the area under f(x^2) between a couple of points if this is the correct interpretation I am not seeing how it will help you find your constants. I think you are closer to your mark in the subject line and may be looking for a "linear-least-squares" fit. (or if you are less lucky a "non-linear-least-squares" fit) disclaimer: I have only been through standard undergraduate classes and have not used this stuff in several years (winter 99, and never in the real world). if I were you I would try to view what I needed as 2 problems each with 3 unknowns since your x and y equations seem independent and the amount of work will normally square with the number of unknowns, that is 2*3^2 < 6^2 looking at just the x half, a few gps samples you could have in something similar to your notation 1*a + lat1*b + lon1*c = x1 1*a + lat2*b + lon2*c = x2 1*a + lat3*b + lon3*c = x3 or in matrix form | 1 lat1 lon1 | |a| |x1| | 1 lat2 lon2 | * |b| = |x2| | 1 lat3 lon3 | |c| |x3| in a more standard notation this is typically represented as A x = b the linear-least-squares solution 'x' to Ax = b is x = (A_t * A)^-1 * A_t * b where A_t is the transpose of A (transpose means the indices for a row and column are swapped or that a matrix is reflected about its diagonal) and the (A_t A)^-1 means the inverse of the matrix that results from multiplying the transpose of A by A (order matters) more specifically it is the matrix you would have to multiply (A_t * A) by to get the "Identity matrix" where the identity matrix is a square matrix with 1's from the top-left to lower-right and 0's elsewhere a matrix must be square to have an inverse There is no guarantee a square matrix has an inverse. I am not going into inverting a matrix because there are alternatives. all in all it is pretty expensive to go this way but at least it is straight forward if you want to plow through it. If you are going to do allot of these there is a shortcut which stems from the rearrangement A_t A x = A_t b where the two sides of the equation can be expressed as Augmented Matrices going back to something that resembles your notation translated to matrix form we get || 1 1 1| |1 lat1 lon1| : a| | 1 1 1 : x1| ||lat1 lat2 lat3| * |1 lat2 lon2| : b| = |lat1 lat2 lat3 : x2| ||lon1 lon2 lon3| |1 lat3 lon3| : c| |lon1 lon2 lon3 : x3| where the [a b c] & [x1 x2 x3] vectors are augmenting the square matrices these augmented matrices (representing your simultaneous equations) may be solved by adding multiples of the rows in a matrix to one another. to end up with 1s, 0s, and solutions (if there are any) there are methods by Gauss, and an improvement called "Gauss-Jordan" which have been around for hundreds of years. in the end you end up with identity matrices augmented with the solutions |1 0 0 : i*a| |1 0 0 : p| |0 1 0 : j*b| = |0 1 0 : q| |0 0 1 : k*c] |0 0 1 : r] so a = p/i b = q/j c = r/k there is a function "naive-gauss" in code I wrote to double check the longhand results required in my homework. (the professor would not even allow a four function calc in the class) and was not finished beyond the point where it served my immediate purpose http://darkwing.uoregon.edu/~tomc/core/math/matrix.r that should work to reduce your augmented equations. there is another "gauss-jordan" which would be better if it did not have a bug I never got around to tracking down note: the "comp" function in the file was derived from Eric Long's compare best of luck On Sat, 9 Nov 2002, [iso-8859-1] jose wrote:

Hi Jose,

> Can you pls send a link to your web article > "Indexless Linear Algebra Algorithms" ? > > The link at numeric-quest seems broken ! > > Thanks again >

Hi Jan, On Sunday, November 10, 2002, 9:14:54 PM, you wrote: JS> Few things have bitten me though. I am still not used JS> to certain "in-place" behaviour of certain functions, such JS> as 'replace, or worst of all 'reverse, which returns empty JS> block. This bites badly in pipelines of expressions. I don't know why REVERSE returns the tail of the block. You can use HEAD REVERSE in expressions, or patch it like this: use [_reverse] [ _reverse: :reverse reverse: func [ "Reverses a series." value [series! tuple! pair!] /part "Limits to a given length or position." range [integer! series!] ] [ head either part [_reverse/part value range] [_reverse value] ] ] Regards, Gabriele. -- Gabriele Santilli <[g--santilli--tiscalinet--it]> -- REBOL Programmer Amigan -- AGI L'Aquila -- REB: http://web.tiscali.it/rebol/index.r

Hi Gabriele,

> I don't know why REVERSE returns the tail of the block. You can > use HEAD REVERSE in expressions, or patch it like this: > ] > ]

Hi Tom, Thanks for helping me out. I think the general approach to solving this is with Jan's code, specially when you want to acount for 4+ data points. It turns out that with 3 data points I only have 2 sets of independent equations, u: f(x,y) and v: f'(x,y), so I only have to solve one set of 3 equations with 3 unknowns and can get the unknows expressed as a function of u,v . Regards, Jose --- Tom Conlin <[tomc--darkwing--uoregon--edu]> escribió:

> Hi Jose, > Hmmm > > > > c, d, e, f below given three calibration > points