Re: Suggestion for dgrid3d

[<HBB...@t-...>]

Philipp K. Janert wrote:
> I still don't think we are really on the same
> page about this.
> 
> We agree that dgrid3d calculates a weighted
> average:
> 	z = Sum_i w_i(d) z_i / Sum_i w_i(d)
> where
> 	z_i is the i-th raw data point 
> 	w_i(d) is the weight function

Let's make that w(d_i).  It's the distance d that depends on i, not the 
function that turns a distance into a w.

We at least partly disagree at this point.  The actual weight function 
is not w(d_i), but

	W_i := w(d_i) / Sum_j w(d_j)

This is smooth (continuous even at d_i = 0), and it's limited to [0:1] 
as weights should be.

> The current implementation uses
> 	w(d) = 1/d**p 	if d > 0
> 	w(d) = 1		if d == 0

No, it doesn't.  It's not really w(d), but the final W_i that's set to 
1.  This is achieved indirectly, by breaking out of the summation loop, 
such that the result becomes the first input point with d==0.

> The lesser reason is this:
> For d small, we multiply values by a large coefficient,
> form the sum, then divide the large coefficient
> out again. Numerically this is not very nice, because
> we loose significant digits if we multiply with something,
> sum, and then divide it out again.

Multiplication doesn't lose digits at all, and summation only does where 
it has to, irrespective of scale.  The multiply/divide pair has 
effectively no effect on numerical accuracy at all.  In other words, the 
given algorithm is quite unlikely to have any bigger problem with 
round-off than a gaussian filter kernel would.