Incremental Averaging

Dear Forthers,

I need to do incremental averaging. I came up with this:

: (average) ( n An Xn+1 -- n+1 An+1 )
swap dup >R
- >R
1 + dup
R> swap /
R> + ;

: average ( seq -- d )
0 0 ['] (average) 2 reduce swap drop ;

REDUCE is a higher order function that takes a stream/"seq", a set of
accumulators, an XT that updates the accumulators and the arity of the
XT. It reads the stream to its end.

Note that in this Forth / also handles floating point division.

My problem is that (AVERAGE) looks really ugly. Any suggestions for
improvement?

Thanks,
Arnold

On 5/8/12 6:45 PM, Arnold Doray wrote:
> Dear Forthers,
>
> I need to do incremental averaging. I came up with this:
>
> : (average) ( n An Xn+1 -- n+1 An+1 )
> swap dup>R
> ->R
> 1 + dup
> R> swap /
> R> + ;
>
> : average ( seq -- d )
> 0 0 ['] (average) 2 reduce swap drop ;
>
> REDUCE is a higher order function that takes a stream/"seq", a set of
> accumulators, an XT that updates the accumulators and the arity of the
> XT. It reads the stream to its end.
>
> Note that in this Forth / also handles floating point division.
>
> My problem is that (AVERAGE) looks really ugly. Any suggestions for
> improvement?

This is the kind of problem M*/ was made for. AVERAGE should accept the
address of an array and its length. Add everything into a double-length
accumulator, and then do an M*/ with a ratio that includes a scale
factor (you know how many actual digits of accuracy you have) and 'n' as
the divisor.

Something like this: (don't have access to a program to test, sorry)

: AVG ( addr n -- n-avg) DUP >R \ Save copy of n
0. 2SWAP 0 DO \ 0. is initial accumulator
>R R@ I CELLS + @ M+ R> LOOP \ Add Ith value to accumulator
10 R> M*/ ; \ Scale average in tenths.

Cheers,
Elizabeth

--
==================================================
Elizabeth D. Rather (US & Canada) 800-55-FORTH
FORTH Inc. +1 310.999.6784
5959 West Century Blvd. Suite 700
Los Angeles, CA 90045
http://www.forth.com

"Forth-based products and Services for real-time
applications since 1973."
==================================================

On Tue, 08 May 2012 21:42:44 -1000, Elizabeth D. Rather wrote:

>
> This is the kind of problem M*/ was made for. AVERAGE should accept the
> address of an array and its length. Add everything into a double-length
> accumulator, and then do an M*/ with a ratio that includes a scale
> factor (you know how many actual digits of accuracy you have) and 'n' as
> the divisor.
>
> Something like this: (don't have access to a program to test, sorry)
>
> : AVG ( addr n -- n-avg) DUP >R \ Save copy of n
> 0. 2SWAP 0 DO \ 0. is initial accumulator
> >R R@ I CELLS + @ M+ R> LOOP \ Add Ith value to accumulator
> 10 R> M*/ ; \ Scale average in tenths.
>
> Cheers,
> Elizabeth

Thank you Elizabeth, but your solution is the plain-vanilla averging, not
incremental averaging. Incremental averaging is less prone to overflow
problems (when you need to read many numbers from a large file or have to
sum squares , eg for RMS or variance calculations).

But thank you for pointing out the use of M*/ . Very useful.

Cheers,
Arnold

Arnold Doray <invalid@invalid.com> writes:
> REDUCE is a higher order function that takes a stream/"seq", a set of
> accumulators, an XT that updates the accumulators and the arity of the
> XT. It reads the stream to its end.

1. That style doesn't seem very Forthish, but whatever.

2. If you just want the final average, it's probably better to add up
everything to a double-width value, and divide at the end. The method
you posted seems to lose precision at every step.

On May 9, 2:54*am, Paul Rubin <no.em...@nospam.invalid> wrote:
> Arnold Doray <inva...@invalid.com> writes:
> > REDUCE is a higher order function that takes a stream/"seq", a set of
> > accumulators, an XT that updates the accumulators and the arity of the
> > XT. It reads the stream to its end.
>
> 1. That style doesn't seem very Forthish, but whatever.

I do stuff like that all of the time in the novice package (EACH
etc.).

What Forth is this?

Maybe I should provide something like your REDUCE in the novice
package. Right now though, I'm confused about what this code does. Is
a "seq" a sequential ascii file with one number on each line? These
are all single-precision integers, right?

> 2. If you just want the final average, it's probably better to add up
> everything to a double-width value, and divide at the end. *The method
> you posted seems to lose precision at every step.

That is what I would do --- simple and straightforward.

Arnold Doray <invalid@invalid.com> writes:
>Dear Forthers,
>
>I need to do incremental averaging. I came up with this:
>
>: (average) ( n An Xn+1 -- n+1 An+1 )
> swap dup >R
> - >R
> 1 + dup
> R> swap /
> R> + ;
[...]
>My problem is that (AVERAGE) looks really ugly. Any suggestions for
>improvement?

With locals:

: (average) { n1 An Xn+1 -- n+1 An+1 }
Xn+1 An - n 1+ tuck / An + ;

Without:

: (average) ( n An Xn+1 -- n+1 An+1 )
over - rot 1+ dup >r / + r> swap ;

Both untested.

- anton
--
M. Anton Ertl http://www.complang.tuwien.ac.at/anton/home.html
comp.lang.forth FAQs: http://www.complang.tuwien.ac.at/forth/faq/toc.html
New standard: http://www.forth200x.org/forth200x.html
EuroForth 2011: http://www.euroforth.org/ef11/

On Wed, 09 May 2012 01:54:12 -0700, Paul Rubin wrote:

> Arnold Doray <invalid@invalid.com> writes:
>> REDUCE is a higher order function that takes a stream/"seq", a set of
>> accumulators, an XT that updates the accumulators and the arity of the
>> XT. It reads the stream to its end.
>
> 1. That style doesn't seem very Forthish, but whatever.

Yes, not mainstream forthish but Reduce, Map, Zip/ZipWith, etc work much
better in Forth than in languages that do support them like Haskell or
Lisp. They are useful because they obviate the need to have explicit
loops, which means they can handle infinite streams.

For example, to sum the first 10,000 squares divsible by 41 :

: ... ( n1 n2 -- seq ) .... ;
: square dup * ;
: pass? 41 mod 0= ;
: sum 0 ['] + 1 REDUCE ;

1 2 ... ' square MAP ' pass? FILTER 10000 TAKE sum

560341738361350000 ok

How to do this cleanly in mainstream Forth? While you can do this in
Haskell, I feel the Forth version is better. YMMV of course.

>
> 2. If you just want the final average, it's probably better to add up
> everything to a double-width value, and divide at the end. The method
> you posted seems to lose precision at every step.

I am dealing with large geophysical datasets, with a few hundred million
entries. I wanted to see the current average as the calculation is being
performed. The incremental version allows this with just a little
tweaking.

Cheers,
Arnold

Arnold Doray <invalid@invalid.com> wrote:
> Dear Forthers,
>
> I need to do incremental averaging. I came up with this:
>
> : (average) ( n An Xn+1 -- n+1 An+1 )
> swap dup >R
> - >R
> 1 + dup
> R> swap /
> R> + ;
>
> : average ( seq -- d )
> 0 0 ['] (average) 2 reduce swap drop ;
>
> REDUCE is a higher order function that takes a stream/"seq", a set of
> accumulators, an XT that updates the accumulators and the arity of the
> XT. It reads the stream to its end.
>
> Note that in this Forth / also handles floating point division.
>
> My problem is that (AVERAGE) looks really ugly. Any suggestions for
> improvement?

It needs a bit of factoring, and I think you're doing too much on the
stack. It's easy if you make the accumulator a structure in memory
with count and average, like so:

0
field: .count
field: .ave
constant /accum

and then:

\ Increment a counter, return its value
: +counter ( a - n) dup @ 1+ tuck ! ;

\ Incremental average
: (average) ( x accum)
dup >r .ave @ -
r@ .count +counter /
r> .ave +! ;

This will also mean less stack thrashing in the word that calls
(AVERAGE) .

Andrew.

On 5/9/12 12:45 AM, Arnold Doray wrote:
> Dear Forthers,
>
> I need to do incremental averaging. I came up with this:
>
> : (average) ( n An Xn+1 -- n+1 An+1 )
> swap dup>R
> ->R
> 1 + dup
> R> swap /
> R> + ;
>
> : average ( seq -- d )
> 0 0 ['] (average) 2 reduce swap drop ;
>
> REDUCE is a higher order function that takes a stream/"seq", a set of
> accumulators, an XT that updates the accumulators and the arity of the
> XT. It reads the stream to its end.
>
> Note that in this Forth / also handles floating point division.
>
> My problem is that (AVERAGE) looks really ugly. Any suggestions for
> improvement?

Interestingly, I came up with locals solution that almost identical to
Anton's, and a non-locals solution that was identical.

In either case, it might help to provide a readable comment such as:

\ An+1 = An + (Xn+1 - An)/(n+1)

Both the locals solution and the comment should help to make this a very
easy-to-maintain program.

-Doug

On Wed, 09 May 2012 03:46:50 -0700, Hugh Aguilar wrote:
>>
>> 1. That style doesn't seem very Forthish, but whatever.
>
> I do stuff like that all of the time in the novice package (EACH etc.).
>
> What Forth is this?
>

It's a Forth I wrote over Java. I am using it now for real apps at work.
It certainly takes the pain out of coding in Java. At least for the
kind of work I do (dealing with geophysical data), this approach looks
promising.

If you are implementing this in "Straight Forth", note that for the
higher order functions to be useful, they have to be lazy, since they act
on infinite streams, a notion I learnt from Clojure.

Some higher order functions I have now are :

REDUCE ( seq <acc> xt <num-cc> -- v )
MAP ( seq xt -- seq' )
ZIP ( seq<a> seq<b> -- seq3<a,b> )
ZIPWITH ( seq<a> seq<b> XT<op> -- seq<op(a,b)> )
TAKE ( seq n -- seq' )
etc.

I don't have folds yet, because I suspect they will be just too slow over
Java.

> Maybe I should provide something like your REDUCE in the novice package.
> Right now though, I'm confused about what this code does. Is a "seq" a
> sequential ascii file with one number on each line? These are all
> single-precision integers, right?

A "seq" is an abstract datastream which has operators "head" and "tail".
It is like a Lisp list, but potentially infinitely long. An idea I stole
from Clojure.

Cheers,
Arnold

Arnold Doray <invalid@invalid.com> writes:
> Yes, not mainstream forthish but Reduce, Map, Zip/ZipWith, etc work much
> better in Forth than in languages that do support them like Haskell or
> Lisp. They are useful because they obviate the need to have explicit
> loops, which means they can handle infinite streams.

Does your Forth have garbage collection? It's been done, but of course
it's not traditional. Without GC it's difficult to use HOF's with much
flexibility. How do you operate on sequences of heap-allocated
structures?

> For example, to sum the first 10,000 squares divsible by 41 : ...
> 560341738361350000 ok

But that answer is wrong, it's actually the sum of the first 100000 of
those squares. Some cut-and-paste error, or a program bug?

> How to do this cleanly in mainstream Forth? While you can do this in
> Haskell, I feel the Forth version is better. YMMV of course.

s = sum $ take 10000 $ [a2 | a <- [1..], let a2 = a*a, a2`mod`41 == 0]

looks natural to me. But say instead of the first 100000 you just want
the ones up to 100000:

t = takeWhile (<= 100000) [a2 | a <- [1..], let a2 = a*a, a2`mod`41 == 0]

and now in addition to the sum, you also want the count:

print (sum t, length t)

How would you do that in your setup? I guess that asks if the right
thing happens if you DUP a sequence.

> I am dealing with large geophysical datasets, with a few hundred million
> entries. I wanted to see the current average as the calculation is being
> performed. The incremental version allows this with just a little
> tweaking.

Well you should probably use floating point, but unless you're worried
about int overflow, keeping a running count and sum seems easier than
that incremental method.

> A "seq" is an abstract datastream which has operators "head" and "tail".
> It is like a Lisp list, but potentially infinitely long. An idea I stole
> from Clojure.

If you want Clojure why don't you use it? Anyway Clojure probably
got it from Scheme, which got it from even older implementations.

anew -average

DOC
(*
> I need to do incremental averaging. I came up with this:
>
> : (average) ( n An Xn+1 -- n+1 An+1 )
> swap dup >R
> - >R
> 1 + dup
> R> swap /
> R> + ;
>
> : average ( seq -- d )
> 0 0 ['] (average) 2 reduce swap drop ;
*)
ENDDOC

: (av) ( An 'Xn+1 n -- An+{Xn+1-An}/n+1 'Xn+2 )
1+ local n+1 @+ local Xn+1 local addr
Xn+1 over - n+1 / + addr ;

CREATE seq #24 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 ,
2 , 2 , 2 , 2 , 2 , 2 , 2 , 2 ,
3 , 3 , 3 , 3 , 3 , 3 , 3 , 3 ,

: average ( addr -- ave ) 0 swap @+ 0 ?DO I (av) LOOP DROP ;

cr .( seq average ) seq average .

: .count ; : .ave cell+ ;
: +counter ( a - n) dup @ 1+ tuck swap ! ;

\ Incremental average
: (average) ( x accum)
dup >r .ave @ -
r@ .count +counter /
r> .ave +! ;

create accum 0 , 0 ,
: 2average ( addr -- ave )
0. accum 2!
@+ 0 ?do @+ accum (average) loop drop
accum .ave @ ;

cr .( seq 2average ) seq 2average .


: (fav) ( F: An -- An+{Xn+1-An}/n+1 ) ( 'Xn+1 n -- 'Xn+2 )
1+ local n+1 @+ local Xn+1 local addr
Xn+1 S>F FOVER F- n+1 S>F F/ F+ addr ;

: faverage ( addr -- ) ( F: -- fave ) 0e @+ 0 ?DO I (fav) LOOP DROP ;

cr .( seq faverage ) seq faverage f.

\ seq average 1
\ seq 2average 1
\ seq faverage 2.000000 ok

-marcel

On Wed, 09 May 2012 11:02:40 -0700, Paul Rubin wrote:
>
> Does your Forth have garbage collection? It's been done, but of course
> it's not traditional. Without GC it's difficult to use HOF's with much
> flexibility. How do you operate on sequences of heap-allocated
> structures?
>

Yes, I built it on top of Java, so it implicitly uses Java's GC. I'm not
sure what you mean by "operate on sequences of heap-allocated
structures".

>> For example, to sum the first 10,000 squares divsible by 41 : ...
>> 560341738361350000 ok
>
> But that answer is wrong, it's actually the sum of the first 100000 of
> those squares. Some cut-and-paste error, or a program bug?
>

My bad. I ran a version on my PC with the first 100,000 squares. Didn't
feel it was fair for most forths, so I took it down a notch but forgot to
update the answer.

>> How to do this cleanly in mainstream Forth? While you can do this in
>> Haskell, I feel the Forth version is better. YMMV of course.
>
> s = sum $ take 10000 $ [a2 | a <- [1..], let a2 = a*a, a2`mod`41 == 0]
>
> looks natural to me. But say instead of the first 100000 you just want
> the ones up to 100000:
>

Haskell is built for this sort of thing. How would you do it in
"traditional" Forth, which was my question.

> t = takeWhile (<= 100000) [a2 | a <- [1..], let a2 = a*a, a2`mod`41 ==
> 0]
>
> and now in addition to the sum, you also want the count:
>
> print (sum t, length t)
>
> How would you do that in your setup? I guess that asks if the right
> thing happens if you DUP a sequence.
>

Yes, the "right thing" certainly happens when you DUP a sequence. That's
the whole point of using seqs. They are not like Java's iterators.

Using :> as my preferred shorthand for :NONAME here's how you would do it
in Forth:

: squares ( seq -- seq' )
['] square map ;

: sum ( seq -- n )
0 ['] + 1 reduce ;

: length ( seq -- n )
0 ['] 1+ 1 reduce ;

1 2 ... squares :> 1000000 <= ; take-while dup sum . length .

333833500
1000001 ok

How would you do this in "mainstream" forth?

>> A "seq" is an abstract datastream which has operators "head" and
>> "tail".
>> It is like a Lisp list, but potentially infinitely long. An idea I
>> stole from Clojure.
>
> If you want Clojure why don't you use it? Anyway Clojure probably got
> it from Scheme, which got it from even older implementations.

Because I need many other things that Clojure doesn't provide.

I needed a number tower (this Forth handles big integers, reals,
rationals & complex numbers), multitasking, REPL, executing on multiple
remote machines, etc. All of which this Forth supports.

The two big takeaways I got from Clojure was the idea of Seqs and that
lazy HOFs with seqs are a good idea, after listening to a talk by its
designer Rich Hickey (It's somewhere on U-Tube, but I can't find a
link).

In any case, I couldn't bring myself to use Clojure on a daily basis.
It's just too ugly. IMHO. But it does have some nice ideas under the
hood.

Forth is capable of presenting the same ideas much more clearly, I feel.
I think the example above proves it to some extent.

Cheers
Arnold

Arnold Doray <invalid@invalid.com> writes:
> Yes, I built it on top of Java, so it implicitly uses Java's GC. I'm not
> sure what you mean by "operate on sequences of heap-allocated
> structures".

Forth stacks traditionally have unboxed integers that are also used as
pointers, characters, etc. If you want non-statically-allocated
multi-word structures they have to be managed somehow, either manually
or with GC. Manual allocation doesn't mix very well with HOF-oriented
programming style, so you want GC and ok, it sounds like you have it.

> Haskell is built for this sort of thing. How would you do it in
> "traditional" Forth, which was my question.

Old fashioned imperative loops and counters.

> Yes, the "right thing" certainly happens when you DUP a sequence. That's
> the whole point of using seqs. They are not like Java's iterators.

Interesting. You probably know about Factor (factorcode.org). I don't
know if it supports that though. I'd say what you're doing is way
outside of normal Forth practice.

> 333833500
> 1000001 ok
> How would you do this in "mainstream" forth?

Loops and stack juggling or accumulation variables.

> I needed a number tower (this Forth handles big integers, reals,
> rationals & complex numbers), multitasking, REPL, executing on multiple
> remote machines, etc. All of which this Forth supports.

I thought Scheme and Common Lisp have that numeric tower and am
surprised if Clojure doesn't (it has a CL compatibility layer, I
thought). I'm told Clojure has good multitasking. I don't know about
the other stuff but I'd expect you can get to the usual Java stuff,
since that's the whole point of running under the JVM.

> Forth is capable of presenting the same ideas much more clearly, I feel.
> I think the example above proves it to some extent.

I think the example is quite far from the traditional concepts of Forth.
That's neither a good nor a bad thing, of course. No language is the
be-all or end-all. Forth's strengh is extreme simplicity of
implementation and predictable execution semantics including timing. I
don't think it's realistic to do anything high speed and hard-real-time
in any JVM environment, but Forth is used for that all the time.

On Thu, 10 May 2012 00:03:01 -0700, Paul Rubin wrote:

>
>> Yes, the "right thing" certainly happens when you DUP a sequence.
>> That's the whole point of using seqs. They are not like Java's
>> iterators.
>
> Interesting. You probably know about Factor (factorcode.org). I don't
> know if it supports that though. I'd say what you're doing is way
> outside of normal Forth practice.
>

Yes, you're right, because "normal" practice is focussed on a much lower
level of absraction - embedded programming, which is Forth's normal
mainstay.

What I saw immediately is that Forth naturally supports much higher
levels , just like Lisp, but far more effectively.

>> 333833500 1000001 ok How would you do this in "mainstream" forth?
>
> Loops and stack juggling or accumulation variables.
>

Of course, Forth is "Turing Complete". However I'd bet the result would
be much harder to debug and understand.

>> I needed a number tower (this Forth handles big integers, reals,
>> rationals & complex numbers), multitasking, REPL, executing on multiple
>> remote machines, etc. All of which this Forth supports.
>
> I thought Scheme and Common Lisp have that numeric tower and am
> surprised if Clojure doesn't (it has a CL compatibility layer, I
> thought). I'm told Clojure has good multitasking. I don't know about
> the other stuff but I'd expect you can get to the usual Java stuff,
> since that's the whole point of running under the JVM.
>

Part of the reason for doing this is that we have a huge load of legacy
code (in Java), which needs to be accessed by non-programmers. I really
don't think they could use Clojure effectively in their work. I can see
them using Forth. It's very simple, and the core language can be easily
extended.

>> Forth is capable of presenting the same ideas much more clearly, I
>> feel.
>> I think the example above proves it to some extent.
>
> I think the example is quite far from the traditional concepts of Forth.
> That's neither a good nor a bad thing, of course. No language is the
> be-all or end-all. Forth's strengh is extreme simplicity of
> implementation and predictable execution semantics including timing. I
> don't think it's realistic to do anything high speed and hard-real-time
> in any JVM environment, but Forth is used for that all the time.

I'm simply saying -- and I hope, proving with real code -- that Forth is
capable of flourishing outside its original design context, and can be
very relevant to solving modern programming problems.

Here is a practical example of reading and parsing a CSV file containing
energy estimates and getting the RMS errors:

: => ( val [name] -- )
create , ;

"scs.csv" open-text => input

\ NB: :> is :NONAME

:> input @ readline ; seq => data

: column ( <string> n -- <string> )
>R " " split to-tuple R> @@ ;

data :> 2 column ; map => E
data :> 3 column ; map => Ei
data :> 4 column ; map => Ew-deficit

: -- ( seq1 seq2 -- seq3 )
' - zipwith ;

Ei Ew-deficit -- => Ew
E Ei -- => Ei-error
E Ew -- => Ew-error

: square ( x -- x*x )
dup * ;

: ^2 ( seq -- seq' )
['] square map ;

: sum ( seq -- n )
0 ['] + 1 reduce ;

: length ( seq -- n )
0 ['] 1+ 1 reduce ;

: average ( seq -- n )
dup sum swap length / ;

Ei-error ^2 average sqrt => RMS-Ei
Ew-error ^2 average sqrt => RMS-Ew

RMS-Ei ? RMS-Ew ?

Cheers,
Arnold