从R中的一组观测数据中生成频率矢量的方法？

Question

问题

我有一个观测向量和它们发生的年份，为了曲线拟合的目的，我想在更长的时间内建立一个频率向量。对于函数，我可以很容易地做到这一点，但是有更简单的方法还是使用固有的矢量化的方法呢？也许我忘记了一些简单的事情。

可复制示例

数据

Events <- data.frame(c(1991, 1991, 1995, 1999, 2007, 2007, 2010, 2010, 2010, 2014), seq(1100, 2000, 100))
names(Events) <- c("Year", "Loss")
Period <- seq(1990, 2014)

函数

FreqV <- function(Period, Observations){
  n <- length(Period)
  F <- double(n)
  for(i in seq_len(n)) {
    F[i] = sum(Observations == Period[i])
  }
  return(F)
}

预期结果

FreqV(Period, Events$Year)
 [1] 0 2 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 2 0 0 3 0 0 0 1

岗位验收更新

这让我感到困扰，为什么算法的C++版本(见“接受答案”下的注释)要慢得多，我终于意识到，原因是上面FreqV的天真翻译。如果有n个周期和m个事件，则必须进行n*m计算。即使在C++，这也是缓慢的。

Tabulate可能被设置为执行一次传递算法，当我在C++中编写一个简单的一次传递算法时，它在5-8之间比表快5-8倍：

天真的C++码

// [[Rcpp::export]]
std::vector<int> FV_C(std::vector<int> P, std::vector<int> O) {
  int n = P.size();
  std::vector<int> F(n);
  for (int i = 0; i < n; ++i){
    F[i] = std::count(O.begin(), O.end(), P[i]);
  }
  return(F);
}

单通C++码

// [[Rcpp::export]]
std::vector<int> FV_C2(std::vector<int> P, std::vector<int> O) {
  int n = P.size();
  int m = O.size();
  int MinP = *std::min_element(P.begin(), P.end());
  std::vector<int> F(n, 0);
  for (int i = 0; i < m; ++i){
    int offset = O[i] - MinP;
    F[offset] +=  1;
  }
  return(F);
}

速度试验

在i7-2600K上使用Windows764bit，使用OpenBLAS 2.13编译的R-3.1.2，用16 an内存在4.6Ghz上进行测试。

set.seed(1)
vals <- sample(sample(10000, 100), 100000, TRUE)
period <- 1:10000

f1a <- function() tabulate(factor(vals, period), nbins = length(period))
f1b <- function() tabulate((vals-period[1])+1, nbins = length(period))
f2 <- function() unname(table(c(period, vals))-1)

library(microbenchmark)

all.equal(f1a(), f1b(), f2(), FV_C(period, vals), FV_C2(period, vals))
[1] TRUE

microbenchmark(f1a(), f1b(), f2(), FV_C(period, vals), FV_C2(period, vals), times = 100L)

Unit: microseconds
                expr        min          lq       mean     median          uq        max neval
               f1a()  26998.194  27812.6250  29515.375  28167.645  28703.4515  55456.079   100
               f1b()    640.049    712.4235   1291.356    800.136   1522.0890  27814.561   100
                f2()  34228.449  35746.6655  39686.660  36210.395  36768.3900  65295.374   100
  FV_C(period, vals) 647577.794 647927.3040 648729.027 648221.417 648848.5090 659463.813   100
 FV_C2(period, vals)    140.877    147.7270    169.085    158.449    170.3625   1095.738   100

Answer 1

我推荐factor和table或tabulate。

这是tabulate

tabulate(factor(Events$Year, Period))
#  [1] 0 2 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 2 0 0 3 0 0 0 1

做这样的事情可能会更快：

tabulate((Events$Year-Period[1])+1)

对于这两种情况，您可能都应该指定nbins，(nbins = length(Period))，以防"Events$Year“中的最大值小于”时间段“中的最大值。

下面是一个性能比较：

set.seed(1)
vals <- sample(sample(10000, 100), 100000, TRUE)
period <- 1:10000

f1a <- function() tabulate(factor(vals, period), nbins = length(period))
f1b <- function() tabulate((vals-period[1])+1, nbins = length(period))
f2 <- function() unname(table(c(period, vals))-1)

library(microbenchmark)
microbenchmark(f1a(), f1b(), f2())
# Unit: microseconds
#   expr       min        lq      mean    median         uq       max neval
#  f1a() 41784.904 43665.394 46789.753 44278.093  45654.546  95032.59   100
#  f1b()   884.465  1162.254  2261.118  1275.154   2756.922  46641.87   100
#   f2() 54837.666 57615.562 71386.516 58863.272 100893.389 130235.33   100

Answer 2

你可以试试

colSums(Vectorize(function(x) x==Events$Year)(Period))
#[1] 0 2 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 2 0 0 3 0 0 0 1

或

colSums(outer(Events$Year, Period, FUN=function(x,y) x==y))
#[1] 0 2 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 2 0 0 3 0 0 0 1

或者使用data.table

library(data.table)
CJ(Period, Events$Year)[, V3:=V1][, sum(V1==V2), V3]$V1
#[1] 0 2 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 2 0 0 3 0 0 0 1

或者如果是命令的话

 c(0,diff(findInterval(Period,Events$Year)))
 #[1] 0 2 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 2 0 0 3 0 0 0 1

或者使用tabulate和fmatch的组合

 library(fastmatch)
 tabulate(fmatch(Events$Year, Period), nbins=length(Period))
 #[1] 0 2 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 2 0 0 3 0 0 0 1

Answer 3

您可以使用table解决这个问题。

table(c(Period,Events$Year))-1

# 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 
#    0    2    0    0    0    1    0    0    0    1    0    0    0    0    0    0    0    2    0    0 
# 2010 2011 2012 2013 2014 
#    3    0    0    0    1 

要去掉这些名字，请使用：

unname(table(c(Period,Events$Year))-1)
# [1] 0 2 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 2 0 0 3 0 0 0 1