Sugar var() only knows numeric

[eddelbuettel]

This fails to compile:

cppFunction('double myVar(IntegerVector x) { return(var(x)); }') 

[wush978]

Dear Dirk,

I study the related header and have some ideas of implementing a template specialization function sugar::Var<INTSXP,NA,T> var( const VectorBase<INTSXP,NA,T>& t) to resolve this issue.

Directly replace the REALSXP to INTSXP will produce the wrong result because the Minus_Vector_Primitive<RTYPE,NA,T> mm( object, m) ; will trim result to int. We need to rewrite the get method.

The simplest approach is to internally transform the IntegerVector to NumericVector and apply the original var.

How do you think?

[eddelbuettel]

But an internal cast would create a copy and that is not ideal.

I think the best way, sadly, is to add four specialization as I did for mean() (cf #219 and #221). It's one of those things I may get to one day when there is nothing more urgent on the plate.

If you want to play with it, this (like a few of the other open issue tickets) is a nice one to create a contribution.

And thanks also for switching to discussion first rather than pull request first ;-)

[wush978]

I'll take this. And thanks for your hint, I'll study your implementation of mean() as a good start.

Ya, discuss first does help!

[wush978]

Dear all,

I don't know if the test is adaquate or not, so I come here for helps. It would be great if any of you can spend 1 minute to review 41fea1e . It compare the result of Rcpp::var on 4 different types with the stats::var.

If the unit test is all what we need, could I send a PR after making the test passed on major machines and compilers?

[eddelbuettel]

Yes, that is the general idea.

In practice most people only use NumericVector which is how we got by so far. But it would be nice to complement it. Not sure anybody ever uses var(myLogicalVector) ... but hey, nice to have.

[MattPD]

Question: Hasn't this commit per chance introduced a regression corresponding to the switch from formula (1) to formula (2) in the following:
http://davegiles.blogspot.com/2015/06/statistical-calculations-numerical.html ?

I'm referring to the following snippet (although CPLXSXP looks analogous):
jmp75@0359729#diff-bfb36f3ad0e86b5efdcbca14acfdaba1L36

I'm wondering, was this intentional?
// The issue is numerical instability & potential negativity: https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance

[MattPD]

Hm, it seems like it has.
Here's a sample snippet demonstrating the issue (edited for clarity):

Rcpp::cppFunction('double var_rcpp_impl(NumericVector x) { return var(x); }') 

vars = function(x)
{
    xbar = mean(x)
    x2 = sum(x ^ 2)
    n = length(x) # EDIT: missing this was the bug
    var_stable = sum((x - xbar) ^ 2) / (n - 1)
    var_unstable = sum(x2 - n * xbar ^ 2) / (n - 1)
    c(var_stable, var_unstable)
}

report_discrepancies = function(x, data_case)
{
    vars_result = vars(x)
    var_R = var(x)
    var_stable = vars_result[1]
    var_unstable = vars_result[2]
    var_rcpp = var_rcpp_impl(x)
    discrepancy_stable = var_R - var_stable
    discrepancy_unstable = var_R - var_unstable
    discrepancy_rcpp = var_R - var_rcpp
    cat("var_R:", var_R, "\n")
    cat(data_case, "; var_stable   =", var_stable, "\t discrepancy = ", discrepancy_stable, "\n")
    cat(data_case, "; var_unstable =", var_unstable, "\t discrepancy = ", discrepancy_unstable, "\n")
    cat(data_case, "; var_rcpp     =", var_rcpp, "\t discrepancy = ", discrepancy_rcpp, "\n")
    cat("\n")
}

set.seed(100)
n = 1e3
x = rnorm(n)

report_discrepancies(x, "benign-data")
report_discrepancies(x + 1e7, "offset-data")
report_discrepancies(x + 1e10, "unpleasant-data")

Output:

var_R: 1.062112 
benign-data ; var_stable   = 1.062112    discrepancy =  0 
benign-data ; var_unstable = 1.062112    discrepancy =  2.220446e-16 
benign-data ; var_rcpp     = 1.062112    discrepancy =  4.440892e-16 

var_R: 1.062112 
offset-data ; var_stable   = 1.062112    discrepancy =  0 
offset-data ; var_unstable = 1.057057    discrepancy =  0.005054458 
offset-data ; var_rcpp     = 1.041041    discrepancy =  0.02107047 

var_R: 1.062111 
unpleasant-data ; var_stable   = 1.062111    discrepancy =  0 
unpleasant-data ; var_unstable = -16794.01   discrepancy =  16795.07 
unpleasant-data ; var_rcpp     = -83970.05   discrepancy =  83971.11 

Following the aforementioned "discussion first" aspect -- I presume changing back to the previous formula (adjusting for the introduced double type instead of STORAGE) could do the trick (EDIT: after reading Thomas Lumley's comment on the previously mentioned blog post -- perhaps we should just do what R does)?

[eddelbuettel]

I am open for suggestions. When I last looked at that I could not create cases with discrepancy.

On the margin, "do what R does but in C++" is a good policy. We may want to be faster at times, but not generally give up precision.

[MattPD]

Right. Will look into comparison with R further (as in: right now I think restoring the previous formula is a strict improvement, but also want to investigate what -- if anything -- breaks if we stop there).

[eddelbuettel]

I am still a little woolly from all the traveling but already 👍 for re-raising this issue. It needs work.

[wush978]

Dear @MattPD and @eddelbuettel ,

Thanks for your great suggestion of the formula. I am glad to help you fix it. I think we need to put a minimal unit test to address this issue and rewrite the implementation of the var function. Your demo is looked great to be a unit test so I'll start from this ASAP.

Wush