Julia で Hotelling の T2 検定

HypothesisTests に Hotelling の T2 検定がある。
分散・共分散行列が等しいと仮定する場合 EqualCovHotellingT2Test と
分散・共分散行列が等しいと仮定しない場合 UnequalCovHotellingT2Test がある。
後者については自由度の計算式がなかなか見つからなかった。
自由度の計算についても，以下の2つのサイトでは違いがある。UnequalCovHotellingT2Test では2番目のサイトと同じであるが，そのサイトでも UnequalCovHotellingT2Test でも，整数自由度で p 値を計算している。

2つの関数を統合して，後者の場合には実数自由度で p 値を計算するようにした。

# https://ncss-wpengine.netdna-ssl.com/wp-content/themes/ncss/pdf/Procedures/NCSS/Hotellings_Two-Sample_T2.pdf

# https://www.real-statistics.com/multivariate-statistics/hotellings-t-square-statistic/hotellings-t-square-unequal-covariance-matrices/

using HypothesisTests
using Rmath, Rmath, LinearAlgebra

function hotellingt2(y1, y2; equalcov=false)
    diff = vec(mean(y1, dims=1) - mean(y2, dims=1));
    n1, p = size(y1);
    n2, q = size(y2);
    p != q && error("column size error");
    if equalcov == true
        spl = ((n1 - 1)cov(y1) + (n2 - 1)cov(y2)) / (n1 + n2 - 2)
        T2 = n1*n2/(n1 + n2) * diff' ⋅ (spl \ diff)
        df2 = n1 + n2 - p - 1
    else
        s1 = cov(y1) / n1;
        s2 = cov(y2) / n2;
        spl = s1 + s2;
        T2 = diff ⋅ (spl \ diff)
        temp = spl \ diff;
        df2 = inv((temp ⋅ (s1 * temp) / T2)^2 / (n1 -1) +
                  (temp ⋅ (s2 * temp) / T2)^2 / (n2 -1))
    end
    F = T2 * (n1 + n2 - p - 1) / (p * (n1 + n2 - 2))
    (T2=T2, F=F, df1=p, df2=df2, pvalue=pf(F, p, df2, false))
end

計算例

テストデータ

y1 = [  49.38   58.892  70.437
        48.668  60.852  72.623
        49.051  60.215  68.984
        50.422  58.335  71.416
        51.423  59.602  72.24
        51.658  61.463  72.736
        49.397  57.89   68.83
        49.442  57.278  62.493
        49.999  61.715  69.557
        50.56   63.758  70.685];

y2 = [  51.548  63.297  68.119
        51.208  64.178  70.406
        52.4    63.338  73.123
        47.647  62.163  68.498
        52.432  63.89   69.763
        51.184  64.765  68.787
        54.983  64.596  74.207
        53.459  64.842  72.325
        52.152  64.65   70.725
        52.091  65.326  73.548
        51.453  62.721  67.53
        53.29   65.564  72.943
        52.192  64.349  71.185
        49.825  64.015  70.478
        52.322  64.005  74.186
        52.443  64.798  71.55
        51.481  63.005  70.572
        51.393  65.698  73.542
        53.494  63.1    66.363
        53.002  61.7    72.151];

分散・共分散行列が等しいと仮定する場合

a = hotellingt2(y1, y2, equalcov=true);
a.T2            # 67.5807155084607
a.F             # 20.91784051452355
a.df1           # 3
a.df2           # 26
a.pvalue        # 4.177382988560072e-7

obj1 = EqualCovHotellingT2Test(y1, y2);
dof(T::EqualCovHotellingT2Test) = (T.p, T.nx + T.ny - T.p - 1)

obj1.T²         # 67.58071550846068
obj1.F          # 20.917840514523544
dof(obj1)       # (3, 26)
pvalue(obj1)    # 4.1773829885600716e-7

分散・共分散行列が等しいと仮定しない場合

b = hotellingt2(y1, y2);
b.T2            # 53.633707938120274
b.F             # 16.60090959989437
b.df1           # 3
b.df2           # 15.163723777577244
b.pvalue        # 4.6975889203673824e-5

obj2 = UnequalCovHotellingT2Test(y1, y2);
dof(T::UnequalCovHotellingT2Test) = (T.p, T.ν)

obj2.T²         # 53.633707938120274
obj2.F          # 16.60090959989437
dof(obj2)       # (3, 15)
pvalue(obj2)    # 4.95288403821627e-5