クラスタリングを評価する過程で、列の分散計算を実行しようとしています (別のトピック)。そうすることで、ベースの「var」関数を使用してこれを行うための素晴らしいトリックに出くわしました。しかし、私が思っていた通りに動作していることを確認するために、手動で分散を計算した結果と比較しました。手巻きバージョンを実行すると、計算でほぼ 1% の差が生じました。浮動小数点の精度による違いがあることは知っていますが、このサイズの数の 1% はちょっと高いようです。
以下は、2つの計算方法をチェックするために使用したコードです（むしろ、3つの方法、モーメントを介して分散も計算し、「var」計算がこのモーメントのショートカットを使用しているように見えます）。サンプル データがあり、分散計算の違いを確認できます。
これらの計算で 1% の差が大きいと断言するのが妥当かどうか疑問に思っています。
ありがとう、
マット
コードは次のとおりです。
> # Sample data
> # This will be 2-d coordinates with pre-defined clusters (treatments)
> dat = rbind(
+ # 4-point cross cluster centered on (2,2)
+ data.frame(grp=1, x=2, y=3),
+ data.frame(grp=1, x=3, y=2),
+ data.frame(grp=1, x=2, y=1),
+ data.frame(grp=1, x=1, y=2),
+ # 2-point dumbbell centered on (-2,-2)
+ data.frame(grp=2, x=-3, y=-2),
+ data.frame(grp=2, x=-1, y=-2),
+ # 3-point equilateral triangle cenetered on (-2,2)
+ data.frame(grp=3, x=-2, y=3),
+ data.frame(grp=3, x=(-2 - sqrt(3)/2), y=-1.5),
+ data.frame(grp=3, x=(-2 + sqrt(3)/2), y=-1.5)
+ )
>
>
> # Compare var calc to hand calc
> # -----------------------------
>
> # Shortcut using existing 'var' function
> var1 = (nrow(dat) - 1) * apply(dat[,-1], 2, var)
> print(var1)
x y
41.05556 37.72222
>
> # Hand-rolled from definition
> centroid = apply(dat[,-1], 2, mean)
> var2 = apply((dat[,-1] - centroid)^2, 2, sum)
> print(var2)
x y
41.46952 37.79630
>
> # Using raw moments
> # Looks to be same as variance calculation
> var3 = apply(dat[,-1], 2, function(col) length(col) * (mean(col^2) - mean(col)^2))
> print(var3)
x y
41.05556 37.72222
>
> # What is the percent difference?
> calcdiff = (var3 - var2) / var2 * 100
> print(calcdiff)
x y
-0.9982268 -0.1959824