SAS proc meansの重み付けの挙動

SASの%stddiffを使って、解析をしたい。
しかし、Rで完結させたい。そのため、前回翻訳をした。

今回、年末年始を利用して最終解析をしようとしたが、その前に
SASで検証してみた。

げげー。なんだか生データから検証すると一部一致しない。

ところが生データからではなく、他の論文中のSDとmeanからは
ズレない。

気持ちが悪いので、調べた。

どうやら、SAS proc meansの重み付けSDがRと
異なっているのが原因。とりあえずRの出来合いの関数
をいくつか試したがどれも一致せず。

また、SASコードを読み解くのか。泣けてくる。
と思ったら、先人のとても参考になるページを発見。感謝。

データステップ100万回　　　　　　SAS新手一生: meansプロシジャのweightステートメントを使用した時の分散とか標準偏差の話
Base SAS(R) 9.2 Procedures Guide

これがベストとは思えないが、現時点のスタンダードなんだろう
として実装しておく。
これで複数の生データの検証でも少なくとも有効数字8桁（SASの
デフォルト。桁の増やし方分かんない）まで一致した。

Githubのアカウントも作ったので、そのうち完全版の%stddiffの
Rのクローンをアップしたいな。

ってか、SDの算出はsqrt(wt.var())の方が良い気がするんだが、
教えてだれか偉い人。

wt.mean <- function (x, wt) #this is a copy of package "SDMTools"
{
      s = which(is.finite(x * wt))
      wt = wt[s]
      x = x[s]
      return(sum(wt * x)/sum(wt))
}

wt.var <- function (x, wt) #this is a copy of package "SDMTools"
{
      s = which(is.finite(x + wt))
      wt = wt[s]
      x = x[s]
      xbar = wt.mean(x, wt)
      return(sum(wt * (x - xbar)^2) * (sum(wt)/(sum(wt)^2 - sum(wt^2))))
}

wt.sd.sas  <- function (x, wt) #Rewriting the code of SAS proc means
{
      s = which(is.finite(x + wt))
      wt = wt[s]
      x = x[s]
      xbar = wt.mean(x, wt)
      return(sqrt(1 / (length(x) - 1) * sum(wt * (x - xbar)^2)))
}