为R中具有相同符号的每个连续数字范围分配一个值
Assigning a value to each range of consecutive numbers with same sign in R
我正在尝试创建一个数据框，其中有一列包含表示正数和负数运行长度的值，如下所示：
Time V Length
0.5 -2 1.5
1.0 -1 1.5
1.5 0 0.0
2.0 2 1.0
2.5 0 0.0
3.0 1 1.75
3.5 2 1.75
4.0 1 1.75
4.5 -1 0.75
5.0 -3 0.75
Length 列对数值为正或负的时间长度求和。零被赋予 0 因为它们是一个拐点。如果没有零分隔符号变化，则对拐点两侧的值进行平均。
我正在尝试估算这些值花费正数或负数的时间量。我已经用 for 循环尝试过这个，并取得了不同程度的成功，但我想避免循环，因为我正在处理非常大的数据集。
我花了一些时间研究 sign 和 diff，因为它们在 this question about sign changes. I've also looked at 中使用，后者使用 transform 和 aggregate 对连续求和重复值。我觉得我可以将它与 sign and/or diff 结合使用，但我不确定如何将这些总和追溯分配给创建它们的范围或如何处理斑点我取拐点处的平均值。
如有任何建议，我们将不胜感激。这是示例数据集：
dat <- data.frame(Time = seq(0.5, 5, 0.5), V = c(-2, -1, 0, 2, 0, 1, 2, 1, -1, -3))
这有效，至少对于您的测试用例而言。它应该非常有效。它做了一些假设，我会尝试指出大的假设。
首先我们提取向量并在开头贴上 0。我们还将最后的V设置为0。计算将基于0s之间的时间差，因此我们需要以0s开始和结束。您的示例似乎默认在 Time = 0 处假设 V = 0，因此初始值为 0，并且它在最长时间突然停止，因此我们也在那里设置 V = 0：
Time = c(0, dat$Time)
V = c(0, dat$V)
V[length(V)] = 0
为了填补跳过的0，我们使用approx对sign(V)进行线性逼近。它还假定您的采样频率是规则的，因此我们可以通过将频率加倍来获得所有缺失的 0。
ap = approx(Time, sign(V), xout = seq(0, max(Time), by = 0.25))
我们要填写的值是观察到的和近似的 0 之间的持续时间。按照正确的顺序，这些是：
dur = diff(ap$x[ap$y == 0])
最后，我们需要原始数据的索引来填充持续时间。这是这个答案中最骇人听闻的部分，但它似乎有效。也许有人会建议一个很好的简化。
# first use rleid to get the sign groupings
group = data.table::rleid(sign(dat$V))
# then we need to set the groups corresponding to 0 values to 0
# and reduce any group numbers following 0s correspondingly
# lastly we add 1 to everything so that we can stick 0 at the
# front of our durations and assign those to the 0 V values
ind = (group - cumsum(dat$V == 0)) * (dat$V != 0) + 1
# fill it in
dat$Length = c(0, dur)[ind]
dat
# Time V Length
# 1 0.5 -2 1.50
# 2 1.0 -1 1.50
# 3 1.5 0 0.00
# 4 2.0 2 1.00
# 5 2.5 0 0.00
# 6 3.0 1 1.75
# 7 3.5 2 1.75
# 8 4.0 1 1.75
# 9 4.5 -1 0.75
# 10 5.0 -3 0.75
首先找到需要插值的"Time"的索引：连续"V"，正负值之间缺少零；他们的 abs(diff(sign(V)) 等于二。
id <- which(abs(c(0, diff(sign(dat$V)))) == 2)
将相关索引之间的平均值为 "Time" 且对应的 "V" 值为零的行添加到原始数据。还要在 "Time" = 0 和最后一个时间步添加 "V" = 0 的行（根据@Gregor 提到的假设）。按 "Time" 排序。
d2 <- rbind(dat,
data.frame(Time = (dat$Time[id] + dat$Time[id - 1])/2, V = 0),
data.frame(Time = c(0, max(dat$Time)), V = c(0, 0))
)
d2 <- d2[order(d2$Time), ]
计算零时间步之间的时间差，并使用 "zero-group indices" 复制它们。
d2$Length <- diff(d2$Time[d2$V == 0])[cumsum(d2$V == 0)]
为原始数据添加值：
merge(dat, d2)
# Time V Length
# 1 0.5 -2 1.50
# 2 1.0 -1 1.50
# 3 1.5 0 1.00
# 4 2.0 2 1.00
# 5 2.5 0 1.75
# 6 3.0 1 1.75
# 7 3.5 2 1.75
# 8 4.0 1 1.75
# 9 4.5 -1 0.75
# 10 5.0 -3 0.75
将 "Length" 设置为 0，其中 V == 0。
我花了比我愿意承认的时间更长的时间，但这是我的解决方案。
因为你说你想在大型数据集上使用它（因此速度很重要）我使用 Rcpp 编写了一个循环来完成所有检查。为了进行速度比较，我还创建了另一个包含 500,000 data.points 的示例数据集并检查了速度（我试图与其他数据集进行比较，但无法将它们转换为 data.table（否则这将是一个不公平的比较...））。如果提供，我会很乐意更新速度比较！
第 1 部分：我的解决方案
我的解决方案如下所示：
(在length_time.cpp)
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
NumericVector length_time(NumericVector time, NumericVector v) {
double start = 0;
double time_i, v_i;
bool last_positive = v[0] > 0;
bool last_negative = v[0] < 0;
int length_i = time.length();
NumericVector ret_vec(length_i);
for (int i = 0; i < length_i; ++i) {
time_i = time[i];
v_i = v[i];
if (v_i == 0) { // injection
if (i > 0) { // if this is not the beginning, then a regime has ended!
ret_vec[i - 1] = time_i - start;
start = time_i;
}
} else if ((v_i > 0 && last_negative) || (v_i < 0 && last_positive)) {
ret_vec[i - 1] = (time_i + time[i - 1]) / 2 - start;
start = (time_i + time[i - 1]) / 2;
}
last_positive = v_i > 0;
last_negative = v_i < 0;
}
ret_vec[length_i - 1] = time[length_i - 1] - start;
// ret_vec now only has the values for the last observation
// do something like a reverse na_locf...
double tmp_val = ret_vec[length_i - 1];
for (int i = length_i - 1; i >= 0; --i) {
if (v[i] == 0) {
ret_vec[i] = 0;
} else if (ret_vec[i] == 0){
ret_vec[i] = tmp_val;
} else {
tmp_val = ret_vec[i];
}
}
return ret_vec;
}
然后在 R 文件中（即 length_time.R）：
library(Rcpp)
# setwd("...") #to find the .cpp-file
sourceCpp("length_time.cpp")
dat$Length <- length_time(dat$Time, dat$V)
dat
# Time V Length
# 1 0.5 -2 1.50
# 2 1.0 -1 1.50
# 3 1.5 0 0.00
# 4 2.0 2 1.00
# 5 2.5 0 0.00
# 6 3.0 1 1.75
# 7 3.5 2 1.75
# 8 4.0 1 1.75
# 9 4.5 -1 0.75
# 10 5.0 -3 0.75
这似乎适用于示例数据集。
第 2 部分：速度测试
library(data.table)
library(microbenchmark)
n <- 10000
set.seed(1235278)
dt <- data.table(time = seq(from = 0.5, by = 0.5, length.out = n),
v = cumsum(round(rnorm(n, sd = 1))))
dt[, chg := v >= 0 & shift(v, 1, fill = 0) <= 0]
plot(dt$time, dt$v, type = "l")
abline(h = 0)
for (i in dt[chg == T, time]) abline(v = i, lty = 2, col = "red")
这会产生一个包含 985 个观测值（交叉点）的数据集。
使用微基准测试速度结果
microbenchmark(dt[, length := length_time(time, v)])
# Unit: milliseconds
# expr min lq mean median uq max neval
# dt[, `:=`(length, length_time(time, v))] 2.625714 2.7184 3.054021 2.817353 3.077489 5.235689 100
计算 500,000 个观测值大约需要 3 毫秒。
对你有帮助吗？
这是我在 base R 中完成的尝试。
Joseph <- function(df) {
is.wholenumber <- function(x, tol = .Machine$double.eps^0.5) abs(x - round(x)) < tol
v <- df$V
t <- df$Time
sv <- sign(v)
nR <- length(v)
v0 <- which(v==0)
id <- which(abs(c(0, diff(sv))) > 1) ## This line and (t[id] + t[id - 1L])/2 From @Henrik
myZeros <- sort(c(v0*t[1L], (t[id] + t[id - 1L])/2))
lenVals <- diff(c(0,myZeros,t[nR])) ## Actual values that
## will populate the Length column
## remove values that result from repeating zeros from the df$V column
lenVals <- lenVals[lenVals != t[1L] | c(!is.wholenumber(myZeros/t[1L]),F)]
## Below we need to determine how long to replicate
## each of the lenVals above, so we need to find
## the starting place and length of each run...
## rle is a great candidate for both of these
m <- rle(sv)
ml <- m$lengths
cm <- cumsum(ml)
zm <- m$values != 0 ## non-zero values i.e. we won't populate anything here
rl <- m$lengths[zm] ## non-zero run-lengths
st <- cm[zm] - rl + 1L ## starting index
out <- vector(mode='numeric', length = nR)
for (i in 1:length(st)) {out[st[i]:(st[i]+rl[i]-1L)] <- lenVals[i]}
df$Length <- out
df
}
这是给定示例的输出：
Joseph(dat)
Time V Length
1 0.5 -2 1.50
2 1.0 -1 1.50
3 1.5 0 0.00
4 2.0 2 1.00
5 2.5 0 0.00
6 3.0 1 1.75
7 3.5 2 1.75
8 4.0 1 1.75
9 4.5 -1 0.75
10 5.0 -3 0.75
这是一个更大的例子：
set.seed(142)
datBig <- data.frame(Time=seq(0.5,50000,0.5), V=sample(-3:3, 10^5, replace=TRUE))
library(compiler)
library(data.table)
library(microbenchmark)
c.Joseph <- cmpfun(Joseph)
c.Henrik <- cmpfun(Henrik)
c.Gregor <- cmpfun(Gregor)
microbenchmark(c.Joseph(datBig), c.Gregor(datBig), c.Henrik(datBig), David(datBig), times = 10)
Unit: milliseconds
expr min lq mean median uq max neval cld
David(datBig) 2.20602 2.617742 4.35927 2.788686 3.13630 114.0674 10 a
c.Joseph(datBig) 61.91015 62.62090 95.44083 64.43548 93.20945 225.4576 10 b
c.Gregor(datBig) 59.25738 63.32861 126.29857 72.65927 214.35961 229.5022 10 b
c.Henrik(datBig) 1511.82449 1678.65330 1727.14751 1730.24842 1816.42601 1871.4476 10 c
正如@Gregor 指出的那样，目标是找到每次出现的零之间的 x 距离。这可以通过绘图直观地看到（再次，正如@Gregor（许多荣誉顺便说一句）指出的那样）。例如，如果我们绘制 datBig 的前 20 个值，我们将获得：
从这里，我们可以看出使图形为正或负（即不为零（当存在重复的零时会发生这种情况））的 x 距离大约为：
2.0, 1.25, 0.5, 0.75, 2.0, 1.0, 0.75, 0.5
t1 <- c.Joseph(datBig)
t2 <- c.Gregor(datBig)
t3 <- c.Henrik(datBig)
t4 <- David(datBig)
## Correct values according to the plot above (x above a value indicates incorrect value)
## 2.00 2.00 2.00 0.00 1.25 1.25 0.50 0.75 0.00 0.00 2.00 2.00 2.00 0.00 0.00 0.00 1.00 0.00 0.75 0.50
## all correct
t1$Length[1:20]
[1] 2.00 2.00 2.00 0.00 1.25 1.25 0.50 0.75 0.00 0.00 2.00 2.00 2.00 0.00 0.00 0.00 1.00 0.00 0.75 0.50
## mostly correct
t2$Length[1:20] x x x x x
[1] 2.00 2.00 2.00 0.00 1.25 1.25 0.50 0.75 0.00 0.00 0.75 0.75 0.75 0.00 0.00 0.00 0.50 0.00 0.75 0.25
## least correct
t3$Length[1:20] x x x x x x x x x x x x x
[1] 2.00 2.00 2.00 0.50 1.00 1.25 0.75 1.25 0.00 1.75 1.75 0.00 1.50 1.50 0.00 0.00 1.25 1.25 1.25 1.25
## all correct
t4$Length[1:20]
[1] 2.00 2.00 2.00 0.00 1.25 1.25 0.50 0.75 0.00 0.00 2.00 2.00 2.00 0.00 0.00 0.00 1.00 0.00 0.75 0.50
# agreement with David's solution
all.equal(t4$Length, t1$Length)
[1] TRUE
嗯，看来 David 提供的 Rcpp 解决方案不仅准确而且速度极快。
我正在尝试创建一个数据框，其中有一列包含表示正数和负数运行长度的值，如下所示：
Time V Length
0.5 -2 1.5
1.0 -1 1.5
1.5 0 0.0
2.0 2 1.0
2.5 0 0.0
3.0 1 1.75
3.5 2 1.75
4.0 1 1.75
4.5 -1 0.75
5.0 -3 0.75
Length 列对数值为正或负的时间长度求和。零被赋予 0 因为它们是一个拐点。如果没有零分隔符号变化，则对拐点两侧的值进行平均。
我正在尝试估算这些值花费正数或负数的时间量。我已经用 for 循环尝试过这个，并取得了不同程度的成功，但我想避免循环，因为我正在处理非常大的数据集。
我花了一些时间研究 sign 和 diff，因为它们在 this question about sign changes. I've also looked at transform 和 aggregate 对连续求和重复值。我觉得我可以将它与 sign and/or diff 结合使用，但我不确定如何将这些总和追溯分配给创建它们的范围或如何处理斑点我取拐点处的平均值。
如有任何建议，我们将不胜感激。这是示例数据集：
dat <- data.frame(Time = seq(0.5, 5, 0.5), V = c(-2, -1, 0, 2, 0, 1, 2, 1, -1, -3))
这有效，至少对于您的测试用例而言。它应该非常有效。它做了一些假设，我会尝试指出大的假设。
首先我们提取向量并在开头贴上 0。我们还将最后的V设置为0。计算将基于0s之间的时间差，因此我们需要以0s开始和结束。您的示例似乎默认在 Time = 0 处假设 V = 0，因此初始值为 0，并且它在最长时间突然停止，因此我们也在那里设置 V = 0：
Time = c(0, dat$Time)
V = c(0, dat$V)
V[length(V)] = 0
为了填补跳过的0，我们使用approx对sign(V)进行线性逼近。它还假定您的采样频率是规则的，因此我们可以通过将频率加倍来获得所有缺失的 0。
ap = approx(Time, sign(V), xout = seq(0, max(Time), by = 0.25))
我们要填写的值是观察到的和近似的 0 之间的持续时间。按照正确的顺序，这些是：
dur = diff(ap$x[ap$y == 0])
最后，我们需要原始数据的索引来填充持续时间。这是这个答案中最骇人听闻的部分，但它似乎有效。也许有人会建议一个很好的简化。
# first use rleid to get the sign groupings
group = data.table::rleid(sign(dat$V))
# then we need to set the groups corresponding to 0 values to 0
# and reduce any group numbers following 0s correspondingly
# lastly we add 1 to everything so that we can stick 0 at the
# front of our durations and assign those to the 0 V values
ind = (group - cumsum(dat$V == 0)) * (dat$V != 0) + 1
# fill it in
dat$Length = c(0, dur)[ind]
dat
# Time V Length
# 1 0.5 -2 1.50
# 2 1.0 -1 1.50
# 3 1.5 0 0.00
# 4 2.0 2 1.00
# 5 2.5 0 0.00
# 6 3.0 1 1.75
# 7 3.5 2 1.75
# 8 4.0 1 1.75
# 9 4.5 -1 0.75
# 10 5.0 -3 0.75
首先找到需要插值的"Time"的索引：连续"V"，正负值之间缺少零；他们的 abs(diff(sign(V)) 等于二。
id <- which(abs(c(0, diff(sign(dat$V)))) == 2)
将相关索引之间的平均值为 "Time" 且对应的 "V" 值为零的行添加到原始数据。还要在 "Time" = 0 和最后一个时间步添加 "V" = 0 的行（根据@Gregor 提到的假设）。按 "Time" 排序。
d2 <- rbind(dat,
data.frame(Time = (dat$Time[id] + dat$Time[id - 1])/2, V = 0),
data.frame(Time = c(0, max(dat$Time)), V = c(0, 0))
)
d2 <- d2[order(d2$Time), ]
计算零时间步之间的时间差，并使用 "zero-group indices" 复制它们。
d2$Length <- diff(d2$Time[d2$V == 0])[cumsum(d2$V == 0)]
为原始数据添加值：
merge(dat, d2)
# Time V Length
# 1 0.5 -2 1.50
# 2 1.0 -1 1.50
# 3 1.5 0 1.00
# 4 2.0 2 1.00
# 5 2.5 0 1.75
# 6 3.0 1 1.75
# 7 3.5 2 1.75
# 8 4.0 1 1.75
# 9 4.5 -1 0.75
# 10 5.0 -3 0.75
将 "Length" 设置为 0，其中 V == 0。
我花了比我愿意承认的时间更长的时间，但这是我的解决方案。
因为你说你想在大型数据集上使用它（因此速度很重要）我使用 Rcpp 编写了一个循环来完成所有检查。为了进行速度比较，我还创建了另一个包含 500,000 data.points 的示例数据集并检查了速度（我试图与其他数据集进行比较，但无法将它们转换为 data.table（否则这将是一个不公平的比较...））。如果提供，我会很乐意更新速度比较！
第 1 部分：我的解决方案
我的解决方案如下所示：
(在length_time.cpp)
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
NumericVector length_time(NumericVector time, NumericVector v) {
double start = 0;
double time_i, v_i;
bool last_positive = v[0] > 0;
bool last_negative = v[0] < 0;
int length_i = time.length();
NumericVector ret_vec(length_i);
for (int i = 0; i < length_i; ++i) {
time_i = time[i];
v_i = v[i];
if (v_i == 0) { // injection
if (i > 0) { // if this is not the beginning, then a regime has ended!
ret_vec[i - 1] = time_i - start;
start = time_i;
}
} else if ((v_i > 0 && last_negative) || (v_i < 0 && last_positive)) {
ret_vec[i - 1] = (time_i + time[i - 1]) / 2 - start;
start = (time_i + time[i - 1]) / 2;
}
last_positive = v_i > 0;
last_negative = v_i < 0;
}
ret_vec[length_i - 1] = time[length_i - 1] - start;
// ret_vec now only has the values for the last observation
// do something like a reverse na_locf...
double tmp_val = ret_vec[length_i - 1];
for (int i = length_i - 1; i >= 0; --i) {
if (v[i] == 0) {
ret_vec[i] = 0;
} else if (ret_vec[i] == 0){
ret_vec[i] = tmp_val;
} else {
tmp_val = ret_vec[i];
}
}
return ret_vec;
}
然后在 R 文件中（即 length_time.R）：
library(Rcpp)
# setwd("...") #to find the .cpp-file
sourceCpp("length_time.cpp")
dat$Length <- length_time(dat$Time, dat$V)
dat
# Time V Length
# 1 0.5 -2 1.50
# 2 1.0 -1 1.50
# 3 1.5 0 0.00
# 4 2.0 2 1.00
# 5 2.5 0 0.00
# 6 3.0 1 1.75
# 7 3.5 2 1.75
# 8 4.0 1 1.75
# 9 4.5 -1 0.75
# 10 5.0 -3 0.75
这似乎适用于示例数据集。
第 2 部分：速度测试
library(data.table)
library(microbenchmark)
n <- 10000
set.seed(1235278)
dt <- data.table(time = seq(from = 0.5, by = 0.5, length.out = n),
v = cumsum(round(rnorm(n, sd = 1))))
dt[, chg := v >= 0 & shift(v, 1, fill = 0) <= 0]
plot(dt$time, dt$v, type = "l")
abline(h = 0)
for (i in dt[chg == T, time]) abline(v = i, lty = 2, col = "red")
这会产生一个包含 985 个观测值（交叉点）的数据集。
使用微基准测试速度结果
microbenchmark(dt[, length := length_time(time, v)])
# Unit: milliseconds
# expr min lq mean median uq max neval
# dt[, `:=`(length, length_time(time, v))] 2.625714 2.7184 3.054021 2.817353 3.077489 5.235689 100
计算 500,000 个观测值大约需要 3 毫秒。
对你有帮助吗？
这是我在 base R 中完成的尝试。
Joseph <- function(df) {
is.wholenumber <- function(x, tol = .Machine$double.eps^0.5) abs(x - round(x)) < tol
v <- df$V
t <- df$Time
sv <- sign(v)
nR <- length(v)
v0 <- which(v==0)
id <- which(abs(c(0, diff(sv))) > 1) ## This line and (t[id] + t[id - 1L])/2 From @Henrik
myZeros <- sort(c(v0*t[1L], (t[id] + t[id - 1L])/2))
lenVals <- diff(c(0,myZeros,t[nR])) ## Actual values that
## will populate the Length column
## remove values that result from repeating zeros from the df$V column
lenVals <- lenVals[lenVals != t[1L] | c(!is.wholenumber(myZeros/t[1L]),F)]
## Below we need to determine how long to replicate
## each of the lenVals above, so we need to find
## the starting place and length of each run...
## rle is a great candidate for both of these
m <- rle(sv)
ml <- m$lengths
cm <- cumsum(ml)
zm <- m$values != 0 ## non-zero values i.e. we won't populate anything here
rl <- m$lengths[zm] ## non-zero run-lengths
st <- cm[zm] - rl + 1L ## starting index
out <- vector(mode='numeric', length = nR)
for (i in 1:length(st)) {out[st[i]:(st[i]+rl[i]-1L)] <- lenVals[i]}
df$Length <- out
df
}
这是给定示例的输出：
Joseph(dat)
Time V Length
1 0.5 -2 1.50
2 1.0 -1 1.50
3 1.5 0 0.00
4 2.0 2 1.00
5 2.5 0 0.00
6 3.0 1 1.75
7 3.5 2 1.75
8 4.0 1 1.75
9 4.5 -1 0.75
10 5.0 -3 0.75
这是一个更大的例子：
set.seed(142)
datBig <- data.frame(Time=seq(0.5,50000,0.5), V=sample(-3:3, 10^5, replace=TRUE))
library(compiler)
library(data.table)
library(microbenchmark)
c.Joseph <- cmpfun(Joseph)
c.Henrik <- cmpfun(Henrik)
c.Gregor <- cmpfun(Gregor)
microbenchmark(c.Joseph(datBig), c.Gregor(datBig), c.Henrik(datBig), David(datBig), times = 10)
Unit: milliseconds
expr min lq mean median uq max neval cld
David(datBig) 2.20602 2.617742 4.35927 2.788686 3.13630 114.0674 10 a
c.Joseph(datBig) 61.91015 62.62090 95.44083 64.43548 93.20945 225.4576 10 b
c.Gregor(datBig) 59.25738 63.32861 126.29857 72.65927 214.35961 229.5022 10 b
c.Henrik(datBig) 1511.82449 1678.65330 1727.14751 1730.24842 1816.42601 1871.4476 10 c
正如@Gregor 指出的那样，目标是找到每次出现的零之间的 x 距离。这可以通过绘图直观地看到（再次，正如@Gregor（许多荣誉顺便说一句）指出的那样）。例如，如果我们绘制 datBig 的前 20 个值，我们将获得：
从这里，我们可以看出使图形为正或负（即不为零（当存在重复的零时会发生这种情况））的 x 距离大约为：
2.0, 1.25, 0.5, 0.75, 2.0, 1.0, 0.75, 0.5
t1 <- c.Joseph(datBig)
t2 <- c.Gregor(datBig)
t3 <- c.Henrik(datBig)
t4 <- David(datBig)
## Correct values according to the plot above (x above a value indicates incorrect value)
## 2.00 2.00 2.00 0.00 1.25 1.25 0.50 0.75 0.00 0.00 2.00 2.00 2.00 0.00 0.00 0.00 1.00 0.00 0.75 0.50
## all correct
t1$Length[1:20]
[1] 2.00 2.00 2.00 0.00 1.25 1.25 0.50 0.75 0.00 0.00 2.00 2.00 2.00 0.00 0.00 0.00 1.00 0.00 0.75 0.50
## mostly correct
t2$Length[1:20] x x x x x
[1] 2.00 2.00 2.00 0.00 1.25 1.25 0.50 0.75 0.00 0.00 0.75 0.75 0.75 0.00 0.00 0.00 0.50 0.00 0.75 0.25
## least correct
t3$Length[1:20] x x x x x x x x x x x x x
[1] 2.00 2.00 2.00 0.50 1.00 1.25 0.75 1.25 0.00 1.75 1.75 0.00 1.50 1.50 0.00 0.00 1.25 1.25 1.25 1.25
## all correct
t4$Length[1:20]
[1] 2.00 2.00 2.00 0.00 1.25 1.25 0.50 0.75 0.00 0.00 2.00 2.00 2.00 0.00 0.00 0.00 1.00 0.00 0.75 0.50
# agreement with David's solution
all.equal(t4$Length, t1$Length)
[1] TRUE
嗯，看来 David 提供的 Rcpp 解决方案不仅准确而且速度极快。