rcpp:在移动窗口计算中删除NAs

Question

我的想法是在一个移动窗口(2乘2)中计算几个统计数据。例如，下面的代码计算移动窗口中的平均值。当输入数据没有NA值时，它工作得很好，但是当NAs在dataset中时，它会给出糟糕的结果(NAs被视为最低的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

这里发生了两件事：

1. 矩阵输入类型arma::imat是一个有符号的int，但是NA和NaN只存在于float或double类型中。本质上，int不能通过设计拥有NA或NaN占位符。因此，发生的转换是放到INT_MIN上。
2. 需要为NA s在C++中子集int或int值。

因此，前进的方法是bale来检测这个INT_MIN值并将其从矩阵中删除。实现这一目标的一种方法是使用find()来识别与INT_MIN和.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