I don't believe this is Kahan summation (which I looked at while trying 
to track this down). Kahan summation seems to work by attempting to 
capture the error in each accumulation to the sum, then adding it back 
to the next element before it is accumulated to the sum.

I believe the R algorithm works as follows.

The first standard calculation of the mean is effectively an estimate of 
the algebraic mean, due to floating point errors (which gets worse the 
further the sum gets away from the elements being accumulated).

The second pass sums the differences of the elements from the estimated 
mean. There should be no net difference as the values on either side of 
the mean should balance out, but we have floating point error. The 
differences from the mean still have the potential for error, but these 
should be smaller than the worst potential difference between a element 
and the accumulating sum (at least the estimated mean lives somewhere 
within the range of values, while the summation may escape it).  
Dividing by N gives you the average difference from the mean, which you 
then use to nudge your initial estimate closer to the true mean. You 
could repeat this to get closer and closer, but at some point the 
floating point error in calculating the average difference from the mean 
will defeat you. I guess one pass is close enough.

All this was explained to me by my wife. Why didn't I think to ask her 
sooner? Choose wisely.





On 2013-07-26 8:58 AM, Ravi Varadhan wrote:
> This uses the idea of Kahan's summation, if I am not mistaken.
>
> http://en.wikipedia.org/wiki/Kahan_summation_algorithm
>
> Ravi
> ________________________________________
> From: r-devel-bounces@r-project.org [r-devel-bounces@r-project.org] on behalf of Joshua Ulrich [josh.m.ulrich@gmail.com]
> Sent: Friday, July 26, 2013 7:12 AM
> To: Zach Harrington
> Cc: r-devel@r-project.org List
> Subject: Re: [Rd] What algorithm is R using to calculate mean?
>
> This was also asked on StackOverflow:
> http://stackoverflow.com/q/17866149/271616.  Here is the answer I
> posted:
>
> This appears to be the updating method of West, 1979 [1] and it was
> implemented in R-2.3.0 in response to PR#1228 [2].
>
> I'm not positive this is the correct algorithm, since it was suggested
> by Martin Maechler but implemented by Brian Ripley. I couldn't find a
> reference in the source code or version control logs that listed the
> actual algorithm used. It was implemented in cov.c in revision 37389
> and in summary.c in revision 37393.
>
> [1] http://dl.acm.org/citation.cfm?doid=359146.359153
> [2] https://bugs.r-project.org/bugzilla3/show_bug.cgi?id=1228
>
> Best,
> --
> Joshua Ulrich  |  about.me/joshuaulrich
> FOSS Trading  |  www.fosstrading.com
>
>
> On Thu, Jul 25, 2013 at 2:44 PM, Zach Harrington
> <zach.harrington@gmail.com> wrote:
>> I am curious to know what algorithm R's mean function uses. Is there some
>> reference to the numerical properties of this algorithm?
>>
>> I found the following C code in summary.c:do_summary():
>> case REALSXP:
>>      PROTECT(ans = allocVector(REALSXP, 1));
>>      for (i = 0; i < n; i++) s += REAL(x)[i];
>>      s /= n;
>>      if(R_FINITE((double)s)) {
>>          for (i = 0; i < n; i++) t += (REAL(x)[i] - s);
>>          s += t/n;
>>      }
>>      REAL(ans)[0] = s;
>>      break;
>>
>> It seems to do a straight up mean:
>> for (i = 0; i < n; i++) s += REAL(x)[i];
>> s /= n;
>>
>> Then it adds what i assume is a numerical correction which seems to be the
>> mean difference from the mean of the data:
>> for (i = 0; i < n; i++) t += (REAL(x)[i] - s);
>> s += t/n;
>>
>> I haven't been able to track this algorithm down anywhere (mean is not a
>> great search term).
>>
>> Any help would be much appreciated,
>>
>> Zach Harrington
>>
>> ______________________________________________
>> R-devel@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-devel
> ______________________________________________
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel


	[[alternative HTML version deleted]]

