Rcpp 精度问题

Question

使用 Rcpp 在 R 中考虑以下 C++ 函数：

cppFunction('long double statZn_cpp(NumericVector dat, double kn) {
  double n = dat.size();
  // Get total sum and sum of squares; this will be the "upper sum"
  // (i.e. the sum above k)
  long double s_upper, s_square_upper;
  // The "lower sums" (i.e. those below k)
  long double s_lower, s_square_lower;
  // Get lower sums
  // Go to kn - 1 to prevent double-counting in main
  // loop
  for (int i = 0; i < kn - 1; ++i) {
    s_lower += dat[i];
    s_square_lower += dat[i] * dat[i];
  }
  // Get upper sum
  for (int i = kn - 1; i < n; ++i) {
    s_upper += dat[i];
    s_square_upper += dat[i] * dat[i];
  }
  // The maximum, which will be returned
  long double M = 0;
  // A candidate for the new maximum, used in a loop
  long double M_candidate;

  // Compute the test statistic
  for (int k = kn; k <= (n - kn); ++k) {
    // Update s and s_square for both lower and upper
    s_lower += dat[k-1];
    s_square_lower += dat[k-1] * dat[k-1];
    s_upper -= dat[k-1];
    s_square_upper -= dat[k-1] * dat[k-1];

    // Get estimate of sd for this k
    long double sdk = sqrt((s_square_lower - pow(s_lower, 2.0) / k +
                      s_square_upper -
                      pow(s_upper, 2.0) / (n - k))/n);
    M_candidate = abs(s_lower / k - s_upper / (n - k)) / sdk;
    // Choose new maximum
    if (M_candidate > M) {
      M = M_candidate;
    }
  }

  return M * sqrt(kn);
}')

尝试命令statZn_cpp(1:20,4)，你会得到6.963106，这是正确答案。缩放应该无关紧要； statZn_cpp(1:20*10,4) 也会得出正确答案 6.963106。但是 statZn_cpp(1:20/10,4) 给出了 6.575959 的错误答案，并且 statZn_cpp(1:20/100,4) 再次给出了明显错误的 0 答案。更重要的是（与我的研究相关，涉及模拟研究），当我尝试 statZn_cpp(rnorm(20),4) 时，答案几乎总是 0，这是错误的。

显然问题与舍入错误有关，但我不知道它们在哪里或如何修复它们（我是 C++ 的新手）。我试图尽可能地扩大精度。有没有办法解决舍入问题？（如果我应该尝试相当于预处理步骤的内容，则允许使用 R 包装函数，但它需要稳健，适用于一般精度水平。）

编辑：这是一些 "equivalent" R 代码：

statZn <- function(dat, kn = function(n) {floor(sqrt(n))}) {
  n = length(dat)
  return(sqrt(kn(n))*max(sapply(
      floor(kn(n)):(n - floor(kn(n))), function(k)
        abs(1/k*sum(dat[1:k]) -
              1/(n-k)*sum(dat[(k+1):n]))/sqrt((sum((dat[1:k] -
               mean(dat[1:k]))^2)+sum((dat[(k+1):n] -
               mean(dat[(k+1):n]))^2))/n))))
}

另外，下面的R代码基本上复制了C++代码应该使用的方法。它能够得出正确答案。

  n = length(dat)
  s_lower = 0
  s_square_lower = 0
  s_upper = 0
  s_square_upper = 0
  for (i in 1:(kn-1)) {
    s_lower = s_lower + dat[i]
    s_square_lower = s_square_lower + dat[i] * dat[i]
  }
  for (i in kn:n) {
    s_upper = s_upper + dat[i]
    s_square_upper = s_square_upper + dat[i] * dat[i]
  }
  M = 0

  for (k in kn:(n-kn)) {
    s_lower = s_lower + dat[k]
    s_square_lower = s_square_lower + dat[k] * dat[k]
    s_upper = s_upper - dat[k]
    s_square_upper = s_square_upper - dat[k] * dat[k]

    sdk = sqrt((s_square_lower - (s_lower)^2/k +
                         s_square_upper -
                         (s_upper)^2/(n-k))/n)
    M_candidate = sqrt(kn) * abs(s_lower / k - s_upper / (n - k)) / sdk

    cat('k', k, '\n',
        "s_lower", s_lower, '\n',
        's_square_lower', s_square_lower, '\n',
        's_upper', s_upper, '\n',
        's_square_upper', s_square_upper, '\n',
        'sdk', sdk, '\n',
        'M_candidate', M_candidate, '\n\n')

    if (M_candidate > M) {
      M = M_candidate
    }
  }

Answer 1

1: 您不应该使用 long double，因为 R 表示 double 类型中的所有数值。使用更精确的类型进行中间计算极不可能提供任何好处，而且更有可能导致平台之间出现奇怪的不一致。

2: 您没有初始化 s_upper、s_square_upper、s_lower 和 s_square_lower。（您实际上是在 R 实现中初始化它们，但您在 C++ 实现中忘记了。）

3: 次要点，但我会用 x*x 替换 pow(x,2.0) 调用。虽然这并不重要。

4: 这就是为我解决的问题：您需要使用包含的命名空间来限定对 C++ 标准库函数的调用。 IOW，std::sqrt() 而不是 sqrt()，std::abs() 而不是 abs()，std::pow() 而不是 pow() 如果你继续使用它.

cppFunction('double statZn_cpp(NumericVector dat, double kn) {
  int n = dat.size();
  double s_upper = 0, s_square_upper = 0; // Get total sum and sum of squares; this will be the "upper sum" (i.e. the sum above k)
  double s_lower = 0, s_square_lower = 0; // The "lower sums" (i.e. those below k)
  for (int i = 0; i < kn - 1; ++i) { s_lower += dat[i]; s_square_lower += dat[i] * dat[i]; } // Get lower sums; Go to kn - 1 to prevent double-counting in main
  for (int i = kn - 1; i < n; ++i) { s_upper += dat[i]; s_square_upper += dat[i] * dat[i]; } // Get upper sum
  double M = 0; // The maximum, which will be returned
  double M_candidate; // A candidate for the new maximum, used in a loop
  // Compute the test statistic
  for (int k = kn; k <= (n - kn); ++k) {
    // Update s and s_square for both lower and upper
    s_lower += dat[k-1];
    s_square_lower += dat[k-1] * dat[k-1];
    s_upper -= dat[k-1];
    s_square_upper -= dat[k-1] * dat[k-1];
    // Get estimate of sd for this k
    double sdk = std::sqrt((s_square_lower - s_lower*s_lower / k + s_square_upper - s_upper*s_upper / (n - k))/n);
    M_candidate = std::abs(s_lower / k - s_upper / (n - k)) / sdk;
    if (M_candidate > M) M = M_candidate; // Choose new maximum
  }
  return std::sqrt(kn) * M;
}');

statZn_cpp(1:20,4); ## you will get 6.963106, which is the correct answer
## [1] 6.963106
statZn_cpp(1:20*10,4); ## Scaling should not matter; will also yield the correct answer of 6.963106
## [1] 6.963106
statZn_cpp(1:20/10,4); ## yields the wrong answer of 6.575959
## [1] 6.963106
statZn_cpp(1:20/100,4); ## again gives you the obviously wrong answer of 0.
## [1] 6.963106
set.seed(1L); statZn_cpp(rnorm(20),4); ## More to the point (and relevant to my research, which involves simulation studies), the answer is almost always 0, which is wrong.
## [1] 1.270117

Rcpp 精度问题

Problems with Rcpp precision

c++

precision

r

rounding

rcpp