What algorithm is R using to calculate mean?
c, numerical-analysis, r
Solution
I'm not sure what algorithm this is, but Martin Maechler mentioned the updating method of West, 1979 in response to PR#1228, which was implemented by Brian Ripley in R-2.3.0. 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.
Problem
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.