Matlabによる初歩の検定

私が使っているMatlabによる検定のプログラムを公開します。 MatlabのStatistics tool boxを前提にしています。私はMatlab 5.3(R11)を使っています。

検定法については「バイオサイエンスの統計学」（南江堂1990）を参考にしました。（なにを見ても検定法はおなじですが）。この本は原理の解説はありませんがそれぞれの検定の使い方を手っ取り早く知りたい人にはお勧めです。また豊富な例は眼力を養うのに役立ちます。

これらのプログラムは無保証ですが、バグを見つけた時はmogami@yugiri.bnf.brain.riken.go.jpにレポートしていただければ感謝します。

t 検定、２標本 t 検定

私は独自のものを使っていますが、これらについてはMatlabのttest, ttest2を使う事をお勧めします。

anov1.m : 一元配置分散分析 (ANOVA, analysis of variance)

function p=anov1(m)
% 1-way ANOVA
% {[group1], [group2], ...}
% each group can have different number of elements
% written by T.Mogami (2003)

k=length(m);
gm=0; %grand mean
gn=0; %total num of samples
for i=1:k
   gm=gm+sum(m{i});
   gn=gn+length(m{i});
end   
gm=gm/gn;
sa=0;
se=0;
for i=1:k
   if length(m{i})~=0
      sa=sa+length(m{i})*(mean(m{i})-gm)^2;
   end
   se=se+ sum((m{i}-gm).^2);      
end
sa=sa/(k-1);
se=se/(gn-k);
p=1-fcdf(sa/se, k-1, gn-k);

Matlabにはanova1関数がありますが、各水準の要素数が等しくなくてはいけないのは、私にとって不都合なので作りました。

contin22.m : ２×２分割表のχ２条独立性の検定 (χ-square test for independence of 2 by 2 contingency table)

function p=contin22(A)
%assume that A is a 2x2 table of contingency
% chi-square test is used.
% It is not accurate, if one of expected frequency is less than 5.

c=sum(A,1);
r=sum(A,2);

N=sum(c);
A
E=r*c /N

chi2=sum(sum((A-E).^2 ./E))
p=1-chi2cdf(chi2, (length(c)-1)*(length(r)-1));

行列として値を与えます。どれかの要素が５を下回るようならこの方法を使わずに、Fisher's exact testを使う必要があります。

gsmooth.m : ガウシアンフィルタ

function t = gsmooth(t, sd)
% gaussian filter with its width sd
% t=0 is assumed outside. filtered in the second dimension of the matrix.
% If the sd is large, this algorithm is not so effective as that use Fourier transformation.
% Written by T.Mogami (2003)

width = 3*ceil(sd);
gw = exp(-((-width):width).^2 /2/sd^2);   %ryougawani 3*sd (~0.3% error)
gw = gw/sum(gw);
t = conv2(t, gw);
t = t(:,(1+width):(size(t,2)+width));

検定ではありませんが、これをほしがる人が多いもので。

kstwo.m : コルモゴロフ=スミルノフテスト

function [prob,d]=kstwo(data1, data2)
%Kolmogolov-Smirnov test
%data1, data2: 1-dimensional array
%return:  p; probability  d; K-S value
%ported from "Numerical recipie in C"
% caution !  This algorithm will not work well if values are degenerate.
j1=1;
j2=1;
fn1=0;
fn2=0;

n1=length(data1);
n2=length(data2);
data1=sort(data1);
data2=sort(data2);
d=0.0;
while(j1<=n1 & j2<=n2)
   d1=data1(j1);
   d2=data2(j2);
   if( d1 <=d2) fn1=j1/n1; j1=j1+1; end
   if( d2 <=d1) fn2=j2/n2; j2=j2+1; end
   dt=abs(fn2-fn1);
   if( dt>d ) d=dt; end
end
en=sqrt(n1*n2/(n1+n2));
prob=probks((en+.12+.11/en)*d);


function sum=probks(alam)
% probability function of Kolmogolov-Smirnov
EPS1=.001;
EPS2=1e-8;
fac=2.0;
sum=.0;
termbf=.0;

a2=-2.0*alam*alam;
for j=1:100
   term=fac*exp(a2*j*j);
   sum=sum+term;
   if (abs(term) <=EPS1*termbf | abs(term)<= EPS2*sum) return; end
   fac=-fac;
   termbf=abs(term);
end
sum=1.0;

Numerical Recipes in C に記述されているアルゴリズムを使っています。

Mann-Whitney検定

Matlabのranksumがそれに当たるのではないかと思いますが確信は持てません。 normalize済みのUを計算したい時は以下の関数を使っています。

function u = MannWhitney(a, b)
% Mann-Whitney U-value, normalized
% assume that a and b are vectors
% 

a=sort(a);
b=sort(b);
s=length(a);
u=0;
for i=s:(-1):1
   u=u+sum(b>a(i))+sum(b==a(i))/2;
end

n1=length(a);
n2=length(b);
mu=n1*n2/2;
sigma=sqrt( n1*n2*(n1+n2+1)/12);
u=(u-mu)/sigma;

Wilcoxon検定

Matlabのsignrankを使いましょう。 ranksumを計算したい時は以下の関数を使っています。

function z = Wilcoxon(a)
% Wilcoxon statistical value T
% assume that a is a vector

[x id] = sort(abs(a));
x = a(id);
n = length(a);
if(sum(x>0) < sum(x<=0))
   T=sum( (x> 0).*(1:n)' )
else
   T=sum( (x<=0).*(1:n)' )
end

m = n*(n+1)/4;
sig = sqrt(n*(n+1)*(2*n+1)/24);
z = (T-m)/sig;

2005年11月19日公開

| Home |
T. Mogami
e-mail: mogami@yugiri.bnf.brain.riken.go.jp