时间戳

首页 > 解决方案 > Rcpp - 如何计算 rowSums 正好为 1 的矩阵

r - Rcpp - 如何计算 rowSums 正好为 1 的矩阵

问题描述

我正在尝试创建一个包含随机数的矩阵，其中rowSums应该正好是 1。

我已经有一个条件来检查是否rowSums不是 1 并尝试更正它。

当我打印出结果时，它看起来是正确的，但是如果我测试所有值是否都是 1，它会给我一些 FALSE 值。

我该如何纠正？

library(Rcpp)

cppFunction('
NumericMatrix imembrandc(int n, int k) {
  NumericMatrix u( n , k );
  IntegerVector sequ = seq(1,100);
  NumericVector sampled;
  for (int i=0; i < k; ++i) {
    sampled = sample(sequ, n);
    u(_,i) = sampled / sum(sampled);
  }

  if (is_true(any(rowSums(u) != 1))) {
    u(_,1) = u(_,1) + (1 - rowSums(u));
  }

  return(u);
}')

当我打印出rowSums结果时，它看起来是正确的：

res = imembrandc(n = 10, k = 5)
rowSums(res)

[1] 1 1 1 1 1 1 1 1 1 1

但是检查它会给出一些错误：

rowSums(res) == 1

[1] 对 对 对 对 对 对 对对对对对对对对对对对对对对对对对对对对对对对对对对对对对对对对对对对对对对对对对对对对对对对对集

标签： rfloating-pointprecisionrcpp

解决方案

生成总和为 1 的随机数的规范方法是从 生成值，将 0 和 1 添加到列表中，并获取排序列表的差。当然，这取决于您想要的随机数分布。这可以用 R 表示为nn - 1[0,1)

set.seed(42)
v <- diff(sort(c(0, runif(5), 1)))
v
#> [1] 0.28613953 0.35560598 0.18870211 0.08435842 0.02226937 0.06292459
sum(v)
#> [1] 1

由reprex 包（v0.2.1）于 2019 年 5 月 24 日创建

在您的 C++ 中：

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericMatrix imembrandc(int n, int k) {
  NumericMatrix u(n, k);
  for (int i = 0; i < n; ++i) {
    NumericVector row = runif(k - 1);
    row.push_back(0.0);
    row.push_back(1.0);
    u(i, _) = diff(row.sort());
  }
  return u;
}

/*** R
set.seed(42)
res = imembrandc(n = 10, k = 5)
rowSums(res)
rowSums(res) == 1
all.equal(rowSums(res),rep(1, nrow(res)))
*/

请注意，当您生成列然后尝试更正rowSum. 输出：

> set.seed(42)

> res = imembrandc(n = 10, k = 5)

> rowSums(res)
 [1] 1 1 1 1 1 1 1 1 1 1

> rowSums(res) == 1
 [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE

> all.equal(rowSums(res),rep(1, nrow(res)))
[1] TRUE

顺便说一句，也all.equal为TRUE您的矩阵提供，因为差异非常小。但我发现最好从一开始就避免这个问题。

推荐阅读