r - 度数分布表から平均偏差と標準偏差を効率的に計算する

Question

次の度数分布表があるとします。

> print(dat)
V1    V2
1  1 11613
2  2  6517
3  3  2442
4  4   687
5  5   159
6  6    29

# V1 = Score
# V2 = Frequency

平均と標準偏差を効率的に計算するにはどうすればよいですか? 収量：SD=0.87 MEAN=1.66。スコアを頻度で複製すると、計算に時間がかかりすぎます。

Answer 1

平均は簡単です。SDは少し注意が必要です（分母にn-1があるため、fastmean（）を再度使用することはできません。

> dat <- data.frame(freq=seq(6),value=runif(6)*100)
> fastmean <- function(dat) {
+   with(dat, sum(freq*value)/sum(freq) )
+ }
> fastmean(dat)
[1] 55.78302
> 
> fastRMSE <- function(dat) {
+   mu <- fastmean(dat)
+   with(dat, sqrt(sum(freq*(value-mu)^2)/(sum(freq)-1) ) )
+ }
> fastRMSE(dat)
[1] 34.9316
> 
> # To test
> expanded <- with(dat, rep(value,freq) )
> mean(expanded)
[1] 55.78302
> sd(expanded)
[1] 34.9316

fastRMSE2回計算することに注意してくださいsum(freq)。これを排除すると、おそらく別のマイナーな速度ブーストが発生します。

ベンチマーク

> microbenchmark(
+   fastmean(dat),
+   mean( with(dat, rep(value,freq) ) )
+   )
Unit: microseconds
                               expr    min      lq median     uq    max
1                     fastmean(dat) 12.433 13.5335 14.776 15.398 23.921
2 mean(with(dat, rep(value, freq))) 21.225 22.3990 22.714 23.406 86.434
> dat <- data.frame(freq=seq(60),value=runif(60)*100)
> 
> dat <- data.frame(freq=seq(60),value=runif(60)*100)
> microbenchmark(
+   fastmean(dat),
+   mean( with(dat, rep(value,freq) ) )
+   )
Unit: microseconds
                               expr    min     lq  median      uq     max
1                     fastmean(dat) 13.177 14.544 15.8860 17.2905  54.983
2 mean(with(dat, rep(value, freq))) 42.610 48.659 49.8615 50.6385 151.053
> dat <- data.frame(freq=seq(600),value=runif(600)*100)
> microbenchmark(
+   fastmean(dat),
+   mean( with(dat, rep(value,freq) ) )
+   )
Unit: microseconds
                               expr      min       lq    median       uq       max
1                     fastmean(dat)   15.706   17.489   25.8825   29.615    79.113
2 mean(with(dat, rep(value, freq))) 1827.146 2283.551 2534.7210 2884.933 26196.923

複製ソリューションは、エントリ数でO（N ^ 2）のように見えます。

このfastmeanソリューションの固定費は12ミリ秒程度で、その後は美しく拡張できます。

より多くのベンチマーク

Comparison with dot product.

dat <- data.frame(freq=seq(600),value=runif(600)*100)
dbaupp <- function(dat) {
  total.count <- sum(dat$freq)
  as.vector(dat$freq %*% dat$value) / total.count
}
microbenchmark(
  fastmean(dat),
  mean( with(dat, rep(value,freq) ) ),
  dbaupp(dat)
)

Unit: microseconds
                               expr     min       lq   median       uq       max
1                       dbaupp(dat)  20.162  21.6875  25.6010  31.3475   104.054
2                     fastmean(dat)  14.680  16.7885  20.7490  25.1765    94.423
3 mean(with(dat, rep(value, freq))) 489.434 503.6310 514.3525 583.2790 30130.302

Answer 2

私は何かが欠けているかもしれませんが、これは非常に迅速に機能するようで、頻度列に数百万を代入することさえできます:

dset <- data.frame(V1=1:6,V2=c(11613,6517,2442,687,159,29))
mean(rep(dset$V1,dset$V2))
#[1] 1.664102
sd(rep(dset$V1,dset$V2))
#[1] 0.8712242

Answer 3

どうですか：

> m = sum(dat$V1 * dat$V2) / sum(dat$V2)
> m
[1] 1.664102
> sigma = sqrt(sum((dat$V1 - m)**2 * dat$V2) / (sum(dat$V2)-1))
> sigma
[1] 0.8712242

ここには複製はありません。

Answer 4

次のコードはレプリケーションを使用せず、R ビルトイン (特にドット積) を可能な限り使用するため、おそらく を使用するソリューションよりも効率的ですsum(V1 * V2)。（編集：これはおそらく間違っています：@ gsk3のソリューションは、私のテストから約1.5〜2倍高速であるようです。）

平均

平均 (または期待値) の定義はsum(n * freq(n)) / total.count、nは「スコア」であり、 は(はちょうど)freq(n)の頻度です。分子の合計は、スコアと度数の内積です。ntotal.countsum(freq(n))

R のドット積は次のとおりです%*%(行列を返しますが、これは基本的にほとんどの目的でベクトルで処理できます)。

> total.count <- sum(dat$V2)
> mean <- dat$V1 %*% dat$V2 / total.count
> mean
         [,1]
[1,] 1.664102

SD

ウィキペディアの記事のこのセクションの最後に数式があり、次のコードに変換されます

> sqrt(dat$V1^2 %*% dat$V2 / total.count - mean^2)
          [,1]
[1,] 0.8712039