Rcpp - Как вычислить матрицу, где rowSums точно 1 - PullRequest
Новая
       29

Rcpp - Как вычислить матрицу, где rowSums точно 1

4 голосов
/ 24 мая 2019

Я пытаюсь создать матрицу со случайными числами, где 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] ИСТИНА ИСТИНА ИСТИНА ИСТИНА ЛОЖЬ ИСТИНА ЛОЖЬ ИСТИНА

1 Ответ

5 голосов
/ 24 мая 2019

Канонический способ для генерации n случайных чисел, сумма которых равна 1, заключается в генерации n - 1 значений из [0,1), добавлении 0 и 1 в список и определении разности отсортированного списка. , Конечно, это зависит от желаемого распределения случайных чисел. Это может быть выражено в R как 

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

Создано в 2019-05-24 с помощью пакета Представление (v0.2.1)

В вашем случае в 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 также для вашей матрицы, поскольку разница действительно небольшая. Но я считаю, что лучше избегать проблемы с самого начала.

...