rcpp：在移动 window 计算中删除 NA

Question

我的想法是在移动window（2乘2）中计算多个统计数据。例如，下面的代码计算移动 window 中的平均值。当输入数据没有 NA 值时它工作得很好，但是当 NA 在数据集中时会给出不好的结果（NA 被视为最低的 int）。你能指导我如何改进它——例如在这些计算中排除 NA 吗？

#include <RcppArmadillo.h>
using namespace Rcpp;

// [[Rcpp::depends(RcppArmadillo)]]

// [[Rcpp::export]]
Rcpp::NumericMatrix get_mw_mean(arma::imat x){
  int num_r = x.n_rows - 1;
  int num_c = x.n_cols - 1;

  arma::dmat result(num_r, num_c);

  for (int i = 0; i < num_r; i++) {
    for (int j = 0; j < num_c; j++) {
      arma::imat sub_x = x.submat(i, j, i + 1, j + 1);
      arma::ivec sub_x_v = vectorise(sub_x);

      arma::vec sub_x_v2 = arma::conv_to<arma::vec>::from(sub_x_v);
      double sub_mean = arma::mean(sub_x_v2);
      result(i, j) = sub_mean;
    }
  }
  return(wrap(result));
}

/*** R
new_c1 = c(1, 86, 98,
           15, 5, 85,
           32, 25, 68)
lg1 = matrix(new_c1, nrow = 3, byrow = TRUE)
get_mw_mean(lg1)
new_c2 = c(NA, 86, 98,
           15, NA, 85,
           32, 25, 68)
lg2 = matrix(new_c2, nrow = 3, byrow = TRUE)
get_mw_mean(lg2)
*/

干杯，记

Answer 1

这里发生了两件事：

矩阵输入类型arma::imat是有符号int，但是NA和NaN 仅出现在 float 或 double 类型中。本质上，int 在设计上不能有 NA 或 NaN 占位符。因此，发生的转换是下降到 INT_MIN。
需要在 C++ 中对 NA 或 NaN 值进行子集化 int。

因此，前进的方向是检测这个 INT_MIN 值并将其从矩阵中删除。实现此目的的一种方法是使用 find() to identify finite elements that do not match INT_MIN and .elem() 来提取已识别的元素。

对于涉及 double 的案例，例如arma::mat/arma::vec/ 等等，考虑使用 find_finite()

已实施

#include <RcppArmadillo.h>

// [[Rcpp::depends(RcppArmadillo)]]

// [[Rcpp::export]]
arma::mat get_mw_mean_na(arma::imat x){
  int num_r = x.n_rows - 1;
  int num_c = x.n_cols - 1;

  Rcpp::Rcout << x <<std::endl;

  arma::dmat result(num_r, num_c);

  for (int i = 0; i < num_r; i++) {
    for (int j = 0; j < num_c; j++) {
      arma::imat sub_x = x.submat(i, j, i + 1, j + 1);
      // Conversion + Search for NA values
      arma::vec sub_x_v2 = arma::conv_to<arma::vec>::from( 
                                        sub_x.elem( find(sub_x != INT_MIN) ) 
      );

      result(i, j) = arma::mean(sub_x_v2);
    }
  }

  return result;
}

输出

new_c1 = c(1, 86, 98,
           15, 5, 85,
           32, 25, 68)
lg1 = matrix(new_c1, nrow = 3, byrow = TRUE)
get_mw_mean_na(lg1)
#        [,1]  [,2]
# [1,] 26.75 68.50
# [2,] 19.25 45.75

new_c2 = c(NA, 86, 98,
           15, NA, 85,
           32, 25, 68)
lg2 = matrix(new_c2, nrow = 3, byrow = TRUE)
get_mw_mean_na(lg2)
#      [,1]     [,2]
# [1,] 50.5 89.66667
# [2,] 24.0 59.33333

rcpp：在移动 window 计算中删除 NA

rcpp: removing NAs in a moving window calculation

r

armadillo

rcpp

已实施

输出