[GRASS5] Re: [GRASS_DE] r.mapcalc round()

[Glynn Clements]

[NB: this thread was last active in mid May 2001]

Glynn Clements wrote:

> > > Does this go away if you compile r.univar with "-ffloat-store"?
> > 
> > Well, r.univar is only a script using r.stats and awk. However, after
> > recompiling r.stats with above flag (adding 
> > EXTRA_CFLAGS=-ffloat-store
> > to the r.stats' Gmakefile) I can't see any difference.
> 
> That wouldn't help; it's the awk script which is doing the relevant
> calculations. Looking at the awk script, I don't think that compiling
> awk with -ffloat-store would help.
> 
> The problem is basically that none of the axioms of real arithmetic
> quite hold true for floating-point arithmetic.
> 
> In this particular case, the variance is:
> 
> 	(SUM[x(i)^2] - SUM[x(i)]^2/N)/N
> 
> If all of the x(i)'s are identical, then the two sides of the
> subtraction should be equal. However, given the error introduced by
> each floating point operation, they will differ slightly; if the RHS
> is slightly larger, then result will be negative, causing the sqrt()
> to fail.
> 
> Whilst there probably is an algorithmic solution (a decent textbook on
> numerical methods should list many approaches for working around the
> problems inherent in floating-point subtraction), I'm not sure whether
> that would be justified, or if it would be better to just include an
> explicit check for a negative variance (indicating a variance so low
> that it's been swamped by rounding error).

FWIW, this subject has just come up on the GMT list, where someone
posted the following URL:

http://web.mit.edu/10.001/Web/Course_Notes/Statistics_Notes/Visualization/node4.html

It describes an algorithm to compute the variance/std.dev. in a single
pass, without the risk of a negative variance due to rounding error. 

-- 
Glynn Clements <glynn.clements at virgin.net>