Hi. This patch reduces the calculation of gsl_histogram_sigma from O(2n) to O(n) and resolves the fixme.
Kind regards, Johannes -- Emails können geändert, gefälscht und eingesehen werden. Signiere oder verschüssele deine Mails mit GPG. http://web.student.tuwien.ac.at/~e0625457/pgp.html
>From 777fab8add4cc83df9c5148b84cad0fd28160e66 Mon Sep 17 00:00:00 2001 From: Johannes Buchner <[email protected]> Date: Sun, 8 Nov 2009 23:41:15 +1300 Subject: [PATCH] a one-pass run algorithm for sigma --- histogram/stat.c | 33 +++++---------------------------- 1 files changed, 5 insertions(+), 28 deletions(-) diff --git a/histogram/stat.c b/histogram/stat.c index 013887d..081ca7d 100644 --- a/histogram/stat.c +++ b/histogram/stat.c @@ -76,15 +76,10 @@ gsl_histogram_sigma (const gsl_histogram * h) const size_t n = h->n; size_t i; - long double wvariance = 0 ; - long double wmean = 0; + long double wsum = 0 ; + long double wsqsum = 0; long double W = 0; - /* FIXME: should use a single pass formula here, as given in - N.J.Higham 'Accuracy and Stability of Numerical Methods', p.12 */ - - /* Compute the mean */ - for (i = 0; i < n; i++) { double xi = (h->range[i + 1] + h->range[i]) / 2; @@ -93,30 +88,12 @@ gsl_histogram_sigma (const gsl_histogram * h) if (wi > 0) { W += wi; - wmean += (xi - wmean) * (wi / W); + wsum += xi * wi; + wsqsum += xi * xi * wi; } } - /* Compute the variance */ - - W = 0.0; - - for (i = 0; i < n; i++) - { - double xi = ((h->range[i + 1]) + (h->range[i])) / 2; - double wi = h->bin[i]; - - if (wi > 0) { - const long double delta = (xi - wmean); - W += wi ; - wvariance += (delta * delta - wvariance) * (wi / W); - } - } - - { - double sigma = sqrt (wvariance) ; - return sigma; - } + return sqrt (1 / W * (wsqsum - wsum * wsum / W)) ; } -- 1.6.4.4
pgpuT16guswXo.pgp
Description: PGP signature
_______________________________________________ Bug-gsl mailing list [email protected] http://lists.gnu.org/mailman/listinfo/bug-gsl
