Merge pull request #320 from MattPD/patch-4

Update var.h

Author: eddelbuettel

ChangeLog

@@ -1,3 +1,11 @@
+2015-07-07  Matt P. Dziubinski  <[email protected]>
+
+        * inst/include/Rcpp/sugar/functions/var.h: Variance -- changed from
+        the unstable formula back to the stable (two-pass) formula, fixed
+        support for complex numbers (formula correction).
+        * inst/unitTests/runit.sugar.var.R: Added tests for complex variance
+        computation, applied simple code refactoring.
+
 2015-07-04  Dirk Eddelbuettel  <[email protected]>
 
         * vignettes/Rcpp.bib: Updated reference to several CRAN packages
@@ -6,7 +14,7 @@
 
         * .clang_format: Added
 
-2015-06-25  Kevin Ushey  <[email protected]>
+2015-06-25 Kevin Ushey  <[email protected]>
 
 	* inst/include/Rcpp/api/meat/Rcpp_eval.h: reset error after Rcpp_eval
 	* inst/unitTests/cpp/Function.cpp: unit tests

inst/NEWS.Rd

@@ -8,10 +8,12 @@
     \itemize{
       \item All internal length variables have been changed from \code{R_len_t}
       to \code{R_xlen_t} to support vectors longer than 2^31-1 elements (via 
-      pull request 303 by Qiang Kou)
-      \item The sugar function \code{sapply} now supports lambda functions
+      pull request 303 by Qiang Kou).
+      \item The sugar function \code{sapply} now supports lambda functions.
+      \item The \code{var} sugar function now uses a more robust two-pass
+      method, supports complex numbers, with new unit tests added.
       \item \code{String} constructors now allow encodings (via pull request 310
-      by Qiang Kou)
+      by Qiang Kou).
     }
     \item Changes in Rcpp Attributes:
     \itemize{

inst/include/Rcpp/sugar/functions/var.h

@@ -34,10 +34,12 @@ class Var : public Lazy< double , Var<RTYPE,NA,T> > {
     Var( const VEC_TYPE& object_ ) : object(object_){}
 
     double get() const{
-        double
-            m = mean(object).get(),
-            ssq = sum( pow(object, 2.0) ).get();
-        return (ssq - m * m * object.size()) / (object.size() - 1);
+        const double average = mean(object).get();
+        const R_xlen_t sample_size = object.size();
+        double sum_squared_deviations = 0.0;
+        for (R_xlen_t i = 0; i != sample_size; ++i)
+            sum_squared_deviations += std::pow(object[i] - average, 2.0);
+        return sum_squared_deviations / (sample_size - 1);
     }
 
 private:
@@ -52,13 +54,14 @@ class Var<CPLXSXP,NA,T> : public Lazy< double , Var<CPLXSXP,NA,T> > {
     Var( const VEC_TYPE& object_ ) : object(object_){}
 
     double get() const{
-        double sq = 0, ssq = 0;
-        for(R_xlen_t i = 0;i < object.size();i++) {
-            Rcomplex z = object[i];
-            sq += z.r;
-            ssq += z.r * z.r;
+        const Rcomplex average = mean(object).get();
+        const R_xlen_t sample_size = object.size();
+        double sum_squared_deviations_magnitudes = 0.0;
+        for (R_xlen_t i = 0; i != sample_size; ++i) {
+            const Rcomplex deviation = object[i] - average;
+            sum_squared_deviations_magnitudes += deviation.r * deviation.r + deviation.i * deviation.i;
         }
-        return (ssq - sq * sq / object.size()) / (object.size() - 1);
+        return sum_squared_deviations_magnitudes / (sample_size - 1);
     }
 
 private:

inst/unitTests/runit.sugar.var.R

@@ -23,14 +23,21 @@
 if (.runThisTest) {
 
     test.Sugar.var <- function() {
-        f1 <- Rcpp::cppFunction('double myVar(NumericVector x) { return(var(x)); }')
-        f2 <- Rcpp::cppFunction('double myVar(IntegerVector x) { return(var(x)); }')
-        f3 <- Rcpp::cppFunction('double myVar(ComplexVector x) { return(var(x)); }')
-        f4 <- Rcpp::cppFunction('double myVar(LogicalVector x) { return(var(x)); }')
-        checkEquals(f1((1:10) * 1.1), var((1:10) * 1.1))
-        checkEquals(f2(1:10), var(1:10))
-        checkEquals(f3(1:10 + (1 + 1i)), var(1:10 + (1 + 1i)))
-        checkEquals(f4(c(T, F, T, F, T)), var(c(T, F, T, F, T)))
+        fNumeric <- Rcpp::cppFunction('double myVar(NumericVector x) { return(var(x)); }')
+        fInteger <- Rcpp::cppFunction('double myVar(IntegerVector x) { return(var(x)); }')
+        fComplex <- Rcpp::cppFunction('double myVar(ComplexVector x) { return(var(x)); }')
+        fLogical <- Rcpp::cppFunction('double myVar(LogicalVector x) { return(var(x)); }')
+        test_data_real <- 1:10
+        checkEquals(fNumeric(test_data_real * 1.1), var(test_data_real * 1.1))
+        checkEquals(fInteger(test_data_real), var(test_data_real))
+        test_data_complex_1 <- complex(real = 5:1, imag = 2:6)
+        test_data_complex_2 <- complex(real = 1:5, imag = 6:10)
+        test_data_complex_1_known_var <- 5
+        test_data_complex_2_known_var <- 5
+        checkEquals(fComplex(test_data_complex_1), test_data_complex_1_known_var)
+        checkEquals(fComplex(test_data_complex_2), test_data_complex_2_known_var)
+        test_data_logical <- c(TRUE, FALSE, TRUE, FALSE, TRUE)
+        checkEquals(fLogical(test_data_logical), var(test_data_logical))
     }
 
 }