计算加权平均数和标准差

Question

我有一个时间序列 x_0 ... x_t。我想计算数据的指数加权方差。也就是说：

V = SUM{w_i*(x_i - x_bar)^2, i=1 to T} where SUM{w_i} = 1 and x_bar=SUM{w_i*x_i}

参考：http://en.wikipedia.org/wiki/Weighted_mean#Weighted_sample_variance

目标是基本上降低时间较早的观测权重。虽然这很容易实现，但我希望尽可能使用内置功能。 有人知道这在R中对应什么吗？

谢谢

- Alex

---

我猜这是一个不完整的规范，你真正想要交付的需要更好地规定如何构建w_i，并且在求和限制方面需要更多细节。 - IRTFM

Answer 1

R提供了加权平均值计算函数。事实上，?weighted.mean展示了以下例子：

 ## GPA from Siegel 1994
 wt <- c(5,  5,  4,  1)/15
 x <- c(3.7,3.3,3.5,2.8)
 xm <- weighted.mean(x, wt)

再来一步：

v <- sum(wt * (x - xm)^2)

- Matthew Lundberg

---

只是一个提示……如果你像我一样迟钝，15是每个重量的总和。然后在这种情况下对重量进行归一化处理。起初我没有理解这一点。 - tharen

你忘记将v除以n或n-1。 - skan

@skan 上述公式表示了该集合的总体方差。请注意，sum(wt) == 1，因此在表达式中乘以wt就相当于除法。 - Matthew Lundberg

它产生的值与SMDtools库中的wt.var(x, wt)函数并不完全相同，也许他们使用了一些修正？ - skan

看起来这个答案目前是最好的。我认为更一致的做法是使用 v <- weighted.mean( (x-xm)^2, wt ) 因为它在权重未被归一化时也能正常工作。 - Michael Lachmann

Answer 2

Hmisc软件包包含您所需的函数。

因此：

x <- c(3.7,3.3,3.5,2.8)

wt <- c(5,  5,  4,  1)/15

xm <- wtd.mean(x, wt)

var <- wtd.var(x, wt)

sd <- sqrt(var)

很遗憾，Hmisc软件包的作者没有包含一个显式的wtd.sd函数。您需要对wtd.var进行平方根运算。

Charles Kangai

- user3770859

---

wtd.mean可以正常工作，但是在你的例子中，wtd.var为“INF”。这是为什么？ - Torvon

@Torvon，这个问题现在已经在Hmisc的开发版本中得到了修复。https://github.com/harrelfe/Hmisc/issues/69 - David J. Harris

sum(wt) 不需要等于1。 - vdesai

Answer 3

当我使用wtd.var()函数时，我也会从Hmisc中得到错误信息。幸运的是，SDMTools具有类似功能，并且甚至可以直接为您计算标准差（SD），而无需取方差的平方根。

library(SDMTools)

x <- c(3.7,3.3,3.5,2.8)
wt <- c(5,  5,  4,  1)/15  ## Note: no actual need to normalize weights to sum to 1, this will be done automatically.

wt.mean(x, wt)
wt.sd(x,wt)

wt.var(x, wt)

- Michael Ohlrogge

---

SDMTools已不再维护。 - Benjamin Ziepert

Answer 4

软件包Hmisc有一个函数wt.var()，正如其他人所指出的。

请注意您需要了解您是否需要频率权重或可靠性权重。在您的情况下，我相信您对可靠性权重感兴趣，因此需要明确设置normwt=TRUE。在这种情况下，您可以以任何格式提供权重（总和为1、总和为N等）。如果您要使用频率权重，就需要小心指定权重的方式。

library(Hmisc)

n <- 3
x <- seq_len(n)
w <- c(0.1, 0.2, 0.6)
w2 <- w / min(w)
w3 <- w / sum(w)

## reliability weights?
wtd.var(x = x, weights = w, normwt=TRUE)
#> [1] 0.95
wtd.var(x = x, weights = w2, normwt=TRUE)
#> [1] 0.95
wtd.var(x = x, weights = w3, normwt=TRUE)
#> [1] 0.95

## frequency weights?
wtd.var(x = x, weights = w)
#> Warning in wtd.var(x = x, weights = w): only one effective observation; variance
#> estimate undefined
#> [1] -4.222222
wtd.var(x = x, weights = w2)
#> [1] 0.5277778
wtd.var(x = x, weights = w3)
#> Warning in wtd.var(x = x, weights = w3): only one effective observation;
#> variance estimate undefined
#> [1] Inf

^{该示例由 Reprex软件包（v0.3.0）于2020年8月26日创建。}

- Matifou

---

频率/可靠性权重之间的差异没有得到足够的关注。+1 - Alex

Answer 5

Hmisc包提供了这个功能：

http://rgm2.lab.nig.ac.jp/RGM2/func.php?rd_id=Hmisc:wtd.stats

- Alex

Answer 6

对于方差和标准差，你必须区分有偏估计和频率、可靠性/抽样权重。

x <- c(1, 4, 5)
wf <- c(1, 2, 3)        # Frequency counts
ws <- c(0.1, 0.2, 0.3)  # Sampling weights

## Weighted mean
mean(rep(x, wf))        # Works only for integer frequencys
#[1] 4
sum(x * wf) / sum(wf)
#[1] 4
sum(x * ws) / sum(ws)
#[1] 4
weighted.mean(x, wf)
#[1] 4
weighted.mean(x, ws)
#[1] 4

## Frequency counts
var(rep(x, wf))                        # Variance
#[1] 2.4
sd(rep(x, wf))                         # Standard deviation
#[1] 1.549193
sw <- sum(wf)
xm <- sum(x * wf) / sw
sum(wf * (x - xm)^2) / (sw - 1)        # Variance
#[1] 2.4
sqrt(sum(wf * (x - xm)^2) / (sw - 1))  # Standard deviation
#[1] 1.549193

## Sampling weights
xm <- weighted.mean(x, ws)
sum(ws *(x-xm)^2)*(sum(ws)/(sum(ws)^2-sum(ws^2)))  # Variance
#[1] 3.272727
cov.wt(matrix(x, ncol=1), ws)$cov                  # Variance
#[1,] 3.272727

## BIASED weighted sample variance
xm <- weighted.mean(x, ws)
sum(ws * (x - xm)^2) / sum(ws)  # Variance
#[1] 2
xm <- weighted.mean(x, wf)
sum(wf * (x - xm)^2) / sum(wf)  # Variance
#[1] 2

## Using Hmisc
Hmisc::wtd.var(x, wf)
#[1] 2.4
Hmisc::wtd.var(x, ws, normwt=TRUE)
#[1] 3.272727

- GKi