计算数据向量的样本统计数据,其存储为频率 table
Compute sample statistics for a data vector with ties which is stored as a frequency table
我正在尝试从具有绑定值的数据向量中获取一些汇总统计信息(均值、方差和分位数)。特别是,它存储在频率分布 table 中:唯一数据值 var
和关系数 frequency
。
我知道我可以使用 rep
函数首先将矢量展开为完整格式:
xx <- rep(mydata$var, mydata$frequency)
然后做标准
mean(xx)
var(xx)
quantile(xx)
但是频率真的很大而且我有很多唯一值,这让程序真的很慢。有没有办法直接从 var
和 frequency
计算这些统计数据?
set.seed(0)
x <- runif(10) ## unique data values
k <- sample.int(5, 10, TRUE) ## frequency
n <- sum(k)
xx <- rep.int(x, k) ## "expanded" data
#################
## sample mean ##
#################
mean(xx) ## using `xx`
#[1] 0.6339458
mu <- c(crossprod(x, k)) / n ## using `x` and `k`
#[1] 0.6339458
#####################
## sample variance ##
#####################
var(xx) * (n - 1) / n ## using `xx`
#[1] 0.06862544
v <- c(crossprod(x ^ 2, k)) / n - mu * mu ## using `x` and `k`
#[1] 0.06862544
计算分位数要复杂得多,但可行。我们需要先了解如何以标准方式计算分位数。
xx <- sort(xx)
pp <- seq(0, 1, length = n)
plot(pp, xx); abline(v = pp, col = 8, lty = 2)
但是,当数据有联系时,我们可以清楚地看到图中有 "runs"(相同值)和 "jumps"(两个值之间)。仅在 "jumps" 上需要线性插值,而在 "runs" 上,分位数只是 运行 值。
以下函数仅使用 x
和 k
查找分位数。出于演示目的,有一个参数 verbose
。如果 TRUE
它将生成包含 "runs"(和 "jumps")信息的绘图和数据框。
find_quantile <- function (x, k, prob = seq(0, 1, length = 5), verbose = FALSE) {
if (is.unsorted(x)) {
ind <- order(x); x <- x[ind]; k <- k[ind]
}
m <- length(x) ## number of unique values
n <- sum(k) ## number of data
d <- 1 / (n - 1) ## break [0, 1] into (n - 1) intervals
## the right and left end of each run
r <- (cumsum(k) - 1) * d
l <- r - (k - 1) * d
if (verbose) {
breaks <- seq(0, 1, d)
plot(r, x, "n", xlab = "prob (p)", ylab = "quantile (xq)", xlim = c(0, 1))
abline(v = breaks, col = 8, lty = 2)
## sketch each run
segments(l, x, r, x, lwd = 3)
## sketch each jump
segments(r[-m], x[-m], l[-1], x[-1], lwd = 3, col = 2)
## sketch `prob`
abline(v = prob, col = 3)
print( data.frame(x, k, l, r) )
}
## initialize the vector of quantiles
xq <- numeric(length(prob))
run <- rbind(l, r)
i <- findInterval(prob, run, rightmost.closed = TRUE)
## odd integers in `i` means that `prob` lies on runs
## quantiles on runs are just run values
on_run <- (i %% 2) != 0
run_id <- (i[on_run] + 1) / 2
xq[on_run] <- x[run_id]
## even integers in `i` means that `prob` lies on jumps
## quantiles on jumps are linear interpolations
on_jump <- !on_run
jump_id <- i[on_jump] / 2
xl <- x[jump_id] ## x-value to the left of the jump
xr <- x[jump_id + 1] ## x-value to the right of the jump
pl <- r[jump_id] ## percentile to the left of the jump
pr <- l[jump_id + 1] ## percentile to the right of the jump
p <- prob[on_jump] ## probability on the jump
## evaluate the line `(pl, xl) -- (pr, xr)` at `p`
xq[on_jump] <- (xr - xl) / (pr - pl) * (p - pl) + xl
xq
}
使用 verbose = TRUE
将函数应用于上面的示例数据得到:
result <- find_quantile(x, k, prob = seq(0, 1, length = 5), TRUE)
# x k l r
#1 0.2016819 4 0.0000000 0.1111111
#2 0.2655087 2 0.1481481 0.1851852
#3 0.3721239 1 0.2222222 0.2222222
#4 0.5728534 4 0.2592593 0.3703704
#5 0.6291140 2 0.4074074 0.4444444
#6 0.6607978 5 0.4814815 0.6296296
#7 0.8966972 1 0.6666667 0.6666667
#8 0.8983897 3 0.7037037 0.7777778
#9 0.9082078 2 0.8148148 0.8518519
#10 0.9446753 4 0.8888889 1.0000000
数据框的每一行都是一个"run"。 x
给出 运行 值,k
是 运行 长度,l
和 r
是 [= 的左右百分位数65=]。图中,"runs"用黑色横线画出。
"jumps" 的信息由一行的 r
、x
值和下一行的 l
、x
值隐含。图中红线画出"jumps"。
垂直绿线表示我们给出的 prob
值。
计算出的分位数是
result
#[1] 0.2016819 0.5226710 0.6607978 0.8983897 0.9446753
与
相同
quantile(xx, names = FALSE)
#[1] 0.2016819 0.5226710 0.6607978 0.8983897 0.9446753
我正在尝试从具有绑定值的数据向量中获取一些汇总统计信息(均值、方差和分位数)。特别是,它存储在频率分布 table 中:唯一数据值 var
和关系数 frequency
。
我知道我可以使用 rep
函数首先将矢量展开为完整格式:
xx <- rep(mydata$var, mydata$frequency)
然后做标准
mean(xx)
var(xx)
quantile(xx)
但是频率真的很大而且我有很多唯一值,这让程序真的很慢。有没有办法直接从 var
和 frequency
计算这些统计数据?
set.seed(0)
x <- runif(10) ## unique data values
k <- sample.int(5, 10, TRUE) ## frequency
n <- sum(k)
xx <- rep.int(x, k) ## "expanded" data
#################
## sample mean ##
#################
mean(xx) ## using `xx`
#[1] 0.6339458
mu <- c(crossprod(x, k)) / n ## using `x` and `k`
#[1] 0.6339458
#####################
## sample variance ##
#####################
var(xx) * (n - 1) / n ## using `xx`
#[1] 0.06862544
v <- c(crossprod(x ^ 2, k)) / n - mu * mu ## using `x` and `k`
#[1] 0.06862544
计算分位数要复杂得多,但可行。我们需要先了解如何以标准方式计算分位数。
xx <- sort(xx)
pp <- seq(0, 1, length = n)
plot(pp, xx); abline(v = pp, col = 8, lty = 2)
以下函数仅使用 x
和 k
查找分位数。出于演示目的,有一个参数 verbose
。如果 TRUE
它将生成包含 "runs"(和 "jumps")信息的绘图和数据框。
find_quantile <- function (x, k, prob = seq(0, 1, length = 5), verbose = FALSE) {
if (is.unsorted(x)) {
ind <- order(x); x <- x[ind]; k <- k[ind]
}
m <- length(x) ## number of unique values
n <- sum(k) ## number of data
d <- 1 / (n - 1) ## break [0, 1] into (n - 1) intervals
## the right and left end of each run
r <- (cumsum(k) - 1) * d
l <- r - (k - 1) * d
if (verbose) {
breaks <- seq(0, 1, d)
plot(r, x, "n", xlab = "prob (p)", ylab = "quantile (xq)", xlim = c(0, 1))
abline(v = breaks, col = 8, lty = 2)
## sketch each run
segments(l, x, r, x, lwd = 3)
## sketch each jump
segments(r[-m], x[-m], l[-1], x[-1], lwd = 3, col = 2)
## sketch `prob`
abline(v = prob, col = 3)
print( data.frame(x, k, l, r) )
}
## initialize the vector of quantiles
xq <- numeric(length(prob))
run <- rbind(l, r)
i <- findInterval(prob, run, rightmost.closed = TRUE)
## odd integers in `i` means that `prob` lies on runs
## quantiles on runs are just run values
on_run <- (i %% 2) != 0
run_id <- (i[on_run] + 1) / 2
xq[on_run] <- x[run_id]
## even integers in `i` means that `prob` lies on jumps
## quantiles on jumps are linear interpolations
on_jump <- !on_run
jump_id <- i[on_jump] / 2
xl <- x[jump_id] ## x-value to the left of the jump
xr <- x[jump_id + 1] ## x-value to the right of the jump
pl <- r[jump_id] ## percentile to the left of the jump
pr <- l[jump_id + 1] ## percentile to the right of the jump
p <- prob[on_jump] ## probability on the jump
## evaluate the line `(pl, xl) -- (pr, xr)` at `p`
xq[on_jump] <- (xr - xl) / (pr - pl) * (p - pl) + xl
xq
}
使用 verbose = TRUE
将函数应用于上面的示例数据得到:
result <- find_quantile(x, k, prob = seq(0, 1, length = 5), TRUE)
# x k l r
#1 0.2016819 4 0.0000000 0.1111111
#2 0.2655087 2 0.1481481 0.1851852
#3 0.3721239 1 0.2222222 0.2222222
#4 0.5728534 4 0.2592593 0.3703704
#5 0.6291140 2 0.4074074 0.4444444
#6 0.6607978 5 0.4814815 0.6296296
#7 0.8966972 1 0.6666667 0.6666667
#8 0.8983897 3 0.7037037 0.7777778
#9 0.9082078 2 0.8148148 0.8518519
#10 0.9446753 4 0.8888889 1.0000000
数据框的每一行都是一个"run"。 x
给出 运行 值,k
是 运行 长度,l
和 r
是 [= 的左右百分位数65=]。图中,"runs"用黑色横线画出。
"jumps" 的信息由一行的 r
、x
值和下一行的 l
、x
值隐含。图中红线画出"jumps"。
垂直绿线表示我们给出的 prob
值。
计算出的分位数是
result
#[1] 0.2016819 0.5226710 0.6607978 0.8983897 0.9446753
与
相同quantile(xx, names = FALSE)
#[1] 0.2016819 0.5226710 0.6607978 0.8983897 0.9446753