Ruby で Mann-Whitney U test をする
Rubyで順位検定のマン- ホイットニーのU検定をするプログラムを書いて見ました．
- 教科書に出てくる平均や分散の近似式で計算しています（配列の長さを n1, n2とすると，平均 n1*n2/2, 分散 sqrt(n1*n2 * (n1+n2-1)/12)
- このため，順位が同じものが非常に多い場合や，2つの配列が両方とも短い場合は，正しい値が求まりません．
- Arrayクラスを拡張する形で実装してあります．
- ブロックには対応していません．
- 同一の順位に対応しようとして，コードがみずらくなったので，キレイにしてくれる人，希望．
使い方
x = [12,11,9,8,8,6,6,4,4,3,3,3,3,2,2,2,1,1,0,0] y = [10,7,7,5,5,4,3,3,2,2,2,2,1,1,1,1,1,1,0,0,0,0,0,0] puts x.mwu(y) # show p-value
- 配列の中は，Float でも大丈夫です．（単に不等号で比較しているだけ）
ソース
- クラスは，Arrayの拡張部分と，z値からp値を求めるメソッド gxpを含むBasicStatFuncの２つです．
- gxpの実装にあたっては，http://aoki2.si.gunma-u.ac.jp/JavaScript/src/gxp.html を移植しました．
#!/usr/bin/env ruby
class BasicStatFunc
def self.gxp(x)
pi2 = 0.398942280401432677940
is = -1.0
y = x.abs
c = y*y
p = 0.0
z = Math.exp(-c*0.5)*pi2
if y < 2.5
20.downto 1 do |i|
p = i*c/(i*2+1+is*p)
is = -is
end
p = 0.5-z*y/(1.0-p)
else
20.downto 1 do |i|
p = i/(y+p)
end
p = z/(y+p)
end
return (x < 0.0) ? 1.0-p : p
end
end
class Array
def mwu(y)
# select short array
return y.mwu(self) if y.size > size
sx = self.sort
sy = y.sort
n = 0
u = 0
yv = sy.shift
sx_size = sx.size
i = 0
while i < sx_size
xv = sx[i]
if yv == xv
st = i
# more than two same values in sx
i += 1 while yv == sx[i]
n = (st + i)/2.0
u += n
# more than two same values in y
while sy.size != 0 && yv == sy[0]
u += n
yv = sy.shift
end
break if sy.size == 0
yv = sy.shift
n = i
elsif yv < xv
u += n
break if sy.size == 0
yv = sy.shift
else
i += 1
n = i
end
end
# when max(y) > max(x)
if i == sx_size
while true
u += n
break if sy.size == 0
yv = sy.shift
end
end
eu = sx_size * y.size
vu = eu * (sx_size + y.size + 1) / 12
eu = eu / 2
z = (u - eu).abs / Math.sqrt(vu)
return BasicStatFunc.gxp(z) * 2 # both sides
end
end
x = [12,11,9,8,8,6,6,4,4,3,3,3,3,2,2,2,1,1,0,0]
y = [10,7,7,5,5,4,3,3,2,2,2,2,1,1,1,1,1,1,0,0,0,0,0,0]
puts "p-value: #{x.mwu(y)}" # 0.034例の値は，バイオサイエンスの統計学から取りました．
バイオサイエンスの統計学―正しく活用するための実践理論
posted with amazlet on 07.05.30
市原 清志
南江堂 (1990/03)
売り上げランキング: 83057
南江堂 (1990/03)
売り上げランキング: 83057
