Mailing List Archive: 49091 messages

# GPS: minimizing sum of squares algorithm

### [1/14] from: jjmmes:y:ahoo:es at: 8-Nov-2002 19:41

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

### [2/14] from: jan:skibinski:sympatico:ca at: 8-Nov-2002 14:22

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:

### [3/14] from: tomc:darkwing:uoregon at: 8-Nov-2002 12:12

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:

### [4/14] from: jjmmes:yaho:o:es at: 9-Nov-2002 16:23

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,

### [5/14] from: jan:skibinski:sympatico:ca at: 9-Nov-2002 11:01

Jose,
> a,b,c,d,e,f are constants. > I basically want to solve a set of 6 equations with
<<quoted lines omitted: 4>>
> Regards > Jose
Ah, then this is quite simple. I thought you used Rebol syntax with equal signs put there by mistake instead of ( :), while in fact you wrote mathematical equations with implied * between 'b and 'lat, etc. Sorry for that. I'll post a solution later today. Jan

### [6/14] 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. }

### [7/14] from: jjmmes::yahoo::es at: 10-Nov-2002 18:11

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,

### [8/14] from: jan:skibinski:sympatico:ca at: 10-Nov-2002 15:14

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! >

### [9/14] from: jjmmes:yaho:o:es at: 10-Nov-2002 21:36

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ó: >

### [11/14] from: jan:skibinski:sympatico:ca at: 10-Nov-2002 16:48

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 >
I know, I broke it. The site is still defunct. But you can retrieve its past contents from the Wayback Machine. Go to www.archive.org, type in a name of a site you are interested in, and voila! Isn't that neat? While in Numeric Quest, select page with Haskell modules. The recommended way is to openit in a separtate window, so you will break out of the stupid frames. All so-called "literate modules", which look like HTML pages, used to be directly executable in Hugs interpreter. Just rename a file giving it a literate extension "lhs" rather than html and you should be able to run it - provided that Wayback Machine did not do too much damage while reformatting the pages - like substituting ">" by >, or by adding java scripts to them. But the Orthogonals.html (from the last archive) seems to be accepted by Hugs in spite of all that reformatting. Enjoy, Jan

### [12/14] from: g:santilli:tiscalinet:it at: 11-Nov-2002 11:31

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

### [13/14] from: jan:skibinski:sympatico:ca at: 11-Nov-2002 11:33

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:
<<quoted lines omitted: 9>>
> ] > ]
Thanks, that's interesting application of 'use! But I'd rather stay away from patching and just get used to the standard behaviour of 'reverse. Jan

### [14/14] from: jjmmes:yah:oo:es at: 11-Nov-2002 21:21

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
<<quoted lines omitted: 153>>
> > > > c, d, e, f below given three calibration > points
=== message truncated ===

Notes
• Quoted lines have been omitted from some messages.
View the message alone to see the lines that have been omitted