[アーティクル]
12/16/2015

November 2015

Volume 30 Number 12

テストの実行 - C# を使用した t 検定

James McCaffrey t 検定は、統計分析の最も基本的な形式の 1 つです。この検定は、標本が 2 セットだけある場合に、この 2 セットの数値の平均が等しくなるかどうかを判断するのが目標です。わかりやすいように例を使って説明します。大きな学区の男子高校生と女子高校生の数学の学力を調べているとします。学力テストには費用も時間もかかるため、すべての生徒にテストを実施することはできません。そこで、標本として 10 人の男子学生と 10 人の女子学生を無作為に選び、数学のテストを実施することにします。この標本の結果に t 検定を実行することで、全男子学生の平均点と全女子学生の平均点が等しくなるかどうかを推測します。

t 検定を実行できるスタンドアロンツールは、Excel を含め、たくさんあります。しかし、スタンドアロンツールを使用する t 検定の機能をソフトウェアシステムに直接統合するのは困難、または不可能です。さらに、著作権などの法的問題が生じる可能性もあります。そこで今回は、外部ライブラリを利用しないで、C# コードだけを使用して、t 検定を実行する方法を説明します。

t 検定の雰囲気を感じ、今回の目標を確認するには、図 1 のデモプログラムを見るのが一番です。最初のデータセットは { 88, 77, 78, 85, 90, 82, 88, 98, 90 } です。これを 10 人の男子学生のテストの点数と考えてください。男子学生の 1 人がなんらかの理由で欠席したため、9 人分の点数だけ記録されています。

図 1 C# を使用した t 検定のデモ

2 つ目のデータセットは { 81, 72, 67, 81, 71, 70, 82, 81 } です。こちらは 10 人の女子学生のテストの点数です。女子学生の 2 人がなんらかの理由で欠席したため、記録が残っているのは 8 人分の点数だけです。最初のデータセットの平均は 86.22 で、2 つ目のデータセットの平均は 75.63 です。2 つの群の平均には約 11 点の差があるため、等しくならないことを示唆しています。しかし、2 つの群 (全男子学生と全女子学生) 全体の平均点が実際には等しかった場合でも、標本データだけを使用することにより、標本の平均点に偶然の差が出る可能性があります。

デモプログラムは、この 2 つの標本データセットを使用して、"t 統計" (t) として値 3.4233 を計算し、"自由度" (df と省略、またはギリシャ文字 ν で表現) として値 14.937 を計算します。次に、求めた t 値と df 値を使用して、確率値 (p-value) として値 0.00379 を計算しています。t 検定の形式はいくつかあります。最も一般的な形式が、スチューデントの t 検定ですが、今回のデモでは、ウェルチの t 検定という改良版を使用します。

p-value は、標本の点数を前提とした場合に、2 つの母集団 (全男子学生と全女子学生) の真の平均点が実際に等しくなる確率です。つまり、計算結果の約 11 点という差が偶然によるものだった確率です。今回の p-value はかなり小さいため、全男子学生の平均点と全女子学生の平均点は等しくならないという結論になります。ほとんど問題では、計算した p-value と比較する境界値を、0.01 にするか 0.05 にするかを自由に定義します。

別の言い方をすると、全男子学生と全女子学生の平均点が等しかった場合、9 人分と 8 人分のデータサイズを持つ 2 つの標本の平均が約 11 点の差になる確率が 0.00379 にすぎないということで、きわめてありそうもない結果といえます。

今回は、少なくとも中級レベルのプログラミングスキルがあることを前提としますが、t 検定の知識は問いません。このデモは C# を使用してコーディングしていますが、コードを Visual Basic .NET や JavaScript などの別の言語にリファクタリングしても大きな問題は起きません。

t 分布について

t 検定は、t 分布に基づきます。t 分布は、正規分布 (ガウス分布または鐘状分布とも呼ばれる) と密接な関係があります。データセットの正規分布の形状は、データの平均と標準偏差によって決まります。標準偏差は、データの散らばり具合 (ばらつき) を測定する 1 つの値です。特殊なケースは、平均 (一般にギリシャ文字 µ) が 0 になるときと、標準偏差 (略語の sd またはギリシャ文字 σ) が 1 になるときです。平均 0 かつ標準偏差 1 の正規分布を、標準正規分布と呼びます。図 2 に、標準正規分布のグラフを示します。

図 2 標準正規分布

図 2 の標準正規分布を定義する方程式を、確率密度関数と呼びます。t 分布は、正規分布とほぼ同じです。t 分布の形状は、"自由度" という 1 つの値によって決まります。図 3 に、自由度 (df) 5 の t 分布を示します。

図 3 の t 分布を定義する方程式には、ギリシャ文字ガンマ (Γ) で示すガンマ関数が関係します。t 検定を実行するには、t 分布の曲線下にある 2 つの等しい領域を計算して合計する必要があります。この 2 つの領域を組み合わせた領域が p-value です。たとえば、図 3 では、t の値が 2.0 の場合、求める組み合わせ領域は、曲線下の負の無限大から -2.0 までの領域と、+2.0 から正の無限大までの領域です。この例の組み合わせ領域 (p-value) は 0.101939 になります。デモプログラムでは、t = 3.4233 となり、この組み合わせ領域が 0.00379 になります。

図 3 t 分布

ここまでは問題ありませんか。では、t 分布の下にある領域をどのように計算するのでしょう。この問題には複数のアプローチがありますが、最も一般的な手法は、標準正規分布の曲線下にある 1 つの関連領域を計算し、その計算結果を使用して p-value を計算することです。たとえば、図 2 では、z (標準偏差の t に相当) が値 -2.0 の場合、負の無限大から -2.0 までの領域を計算し、0.02275 を求めます。次に、正規曲線のこの領域から、t 分布の対応する領域を計算します。

つまり、t 検定を実行するには、t 分布の 2 つの (等しい) 領域を計算し、合計する必要があります。この領域を p-value と呼びます。そのためには、標準正規分布の 1 つの領域を計算後、その結果領域を使用して p-value を求めます。

標準正規分布の領域の計算

標準正規分布曲線下の領域を計算する方法はたくさんあります。これは、コンピューターサイエンスにおける最古の問題の 1 つです。個人的な好みで、いわゆる ACM Algorithm #209 を使用します。Association for Computing Machinery (ACM) は、数値計算と統計計算の基礎アルゴリズムを数多く公開しています。

図 4 に、Algorithm #209 の C# 実装を Gauss 関数として示します。この関数は、負の無限大と正の無限大の間の値 z を受け取り、標準正規分布下にある負の無限大から z までの近接する近似領域を返します。

図 4 標準正規分布下の領域の計算

public static double Gauss(double z)
{
  // input = z-value (-inf to +inf)
  // output = p under Standard Normal curve from -inf to z
  // e.g., if z = 0.0, function returns 0.5000
  // ACM Algorithm #209
  double y; // 209 scratch variable
  double p; // result. called 'z' in 209
  double w; // 209 scratch variable
  if (z == 0.0)
    p = 0.0;
  else
  {
    y = Math.Abs(z) / 2;
    if (y >= 3.0)
    {
      p = 1.0;
    }
    else if (y < 1.0)
    {
      w = y * y;
      p = ((((((((0.000124818987 * w
        - 0.001075204047) * w + 0.005198775019) * w
        - 0.019198292004) * w + 0.059054035642) * w
        - 0.151968751364) * w + 0.319152932694) * w
        - 0.531923007300) * w + 0.797884560593) * y * 2.0;
    }
    else
    {
      y = y - 2.0;
      p = (((((((((((((-0.000045255659 * y
        + 0.000152529290) * y - 0.000019538132) * y
        - 0.000676904986) * y + 0.001390604284) * y
        - 0.000794620820) * y - 0.002034254874) * y
        + 0.006549791214) * y - 0.010557625006) * y
        + 0.011630447319) * y - 0.009279453341) * y
        + 0.005353579108) * y - 0.002141268741) * y
        + 0.000535310849) * y + 0.999936657524;
    }
  }
  if (z > 0.0)
    return (p + 1.0) / 2;
  else
    return (1.0 - p) / 2;
}

図 4 のコードを一目見れば、ゼロから独自の実装を記述するよりも、ACM #209 などの既存のアルゴリズムを使用する方が簡単だと納得するでしょう。ACM #209 の代わりに、Milton Abramowitz および Irene A. Stegun 共著の『Handbook of Mathematical Functions』(Dover Publications、1965 年) に記載されている 7.1.26 方程式を少し編集して使用することもできます。

t 分布の領域の計算

Gauss 関数を実装したら、t 分布の領域を ACM Algorithm #395 を使用して計算します。図 5 に、Algorithm #395 の C# 実装を Student 関数として示します。この関数は t 値と df 値を受け取り、負の無限大から t までと、t から正の無限大までを組み合わせた領域を返します。

図 5 t 分布下の領域の計算

public static double Student(double t, double df)
{
  // for large integer df or double df
  // adapted from ACM algorithm 395
  // returns 2-tail p-value
  double n = df; // to sync with ACM parameter name
  double a, b, y;
  t = t * t;
  y = t / n;
  b = y + 1.0;
  if (y > 1.0E-6) y = Math.Log(b);
  a = n - 0.5;
  b = 48.0 * a * a;
  y = a * y;
  y = (((((-0.4 * y - 3.3) * y - 24.0) * y - 85.5) /
    (0.8 * y * y + 100.0 + b) + y + 3.0) / b + 1.0) *
    Math.Sqrt(y);
  return 2.0 * Gauss(-y); // ACM algorithm 209
}

Algorithm #395 には、2 つの形式があります。1 つは、df パラメーターを整数値として受け取る形式で、もう 1 つは Double 型の値として受け取る形式です。ほとんどの統計計算では自由度を整数値にしますが、ウェルチの t 検定では Double 型の値を使用します。

デモプログラム

デモプログラムを作成するには、Visual Studio を起動して、TTest という名前の新しい C# コンソールアプリケーションを作成します。このデモは特定の .NET バージョンとの依存関係がないので、どのバージョンの Visual Studio でも動作します。テンプレートコードをエディターに読み込んだ後、最上位レベルの System 名前空間への 1 つの参照を除いて、すべての using ステートメントを削除します。次に、ソリューションエクスプローラーウィンドウで、Program.cs ファイルの名前を TTestProgram.cs に変更します。これにより Program クラスの名前が Visual Studio によって自動的に変更されます。

デモプログラムは長すぎてここにすべて掲載することはできませんが、完全なソースコードは本稿付属のコードダウンロードから入手できます。Main メソッドでは、まず、2 つの標本データセットを設定および表示します。

Console.WriteLine("\nBegin Welch's t-test using C# demo\n");
var x = new double[] { 88, 77, 78, 85, 90, 82, 88, 98, 90 };
var y = new double[] { 81, 72, 67, 81, 71, 70, 82, 81 };
Console.WriteLine("\nThe first data set (x) is:\n");
ShowVector(x, 0);
Console.WriteLine("\nThe second data set (y) is:\n");
ShowVector(y, 0);

すべての処理は、TTest というメソッドで実行します。

Console.WriteLine("\nStarting Welch's t-test using C#\n");
TTest(x, y);
Console.WriteLine("\nEnd t-test demo\n");
Console.ReadLine();

TTest メソッドの定義では、まず、各データセットの合計を求めます。

public static void TTest(double[] x, double[] y)
{
  double sumX = 0.0;
  double sumY = 0.0;
  for (int i = 0; i < x.Length; ++i)
    sumX += x[i];
  for (int i = 0; i < y.Length; ++i)
    sumY += y[i];
...

次に、その合計値を使用して、2 つの標本の平均を計算します。

int n1 = x.Length;
int n2 = y.Length;
double meanX = sumX / n1;
double meanY = sumY / n2;

続いて、2 つの平均を使用して、2 つの標本の分散を計算します。

double sumXminusMeanSquared = 0.0; // Calculate variances
double sumYminusMeanSquared = 0.0;
for (int i = 0; i < n1; ++i)
  sumXminusMeanSquared += (x[i] - meanX) * (x[i] - meanX);
for (int i = 0; i < n2; ++i)
  sumYminusMeanSquared += (y[i] - meanY) * (y[i] - meanY);
double varX = sumXminusMeanSquared / (n1 - 1);
double varY = sumYminusMeanSquared / (n2 - 1);

データ集合の分散は標準偏差の平方になります。したがって、標準偏差は分散の平方根です。t 検定はこの分散を使って機能します。ここで t 統計を計算します。

double top = (meanX - meanY);
double bot = Math.Sqrt((varX / n1) + (varY / n2));
double t = top / bot;

言葉で説明すると、まず 2 つの標本の平均値の差を求めます。次に各標本の分散をその標本のサイズで除算したものを合計し、その平方根を求めます。最後に平均値の差を平方根値で除算したものが t 統計です。次に、自由度を計算します。

double num = ((varX / n1) + (varY / n2)) *
  ((varX / n1) + (varY / n2));
double denomLeft = ((varX / n1) * (varX / n1)) / (n1 - 1);
double denomRight = ((varY / n2) * (varY / n2)) / (n2 - 1);
double denom = denomLeft + denomRight;
double df = num / denom;

ウェルチ t 検定の自由度の計算はやや難解で、その方程式はあまり明確ではありません。さいわい、上記の計算に手を加える必要はありません。TTest メソッドは、p-value を計算し、計算した値をすべて表示することによって完結します。

...
  double p = Student(t, df); // Cumulative two-tail density
  Console.WriteLine("mean of x = " + meanX.ToString("F2"));
  Console.WriteLine("mean of y = " + meanY.ToString("F2"));
  Console.WriteLine("t = " + t.ToString("F4"));
  Console.WriteLine("df = " + df.ToString("F3"));
  Console.WriteLine("p-value = " + p.ToString("F5"));
  Explain();
}

Explain というプログラム定義のメソッドは、p-value の解釈方法を説明する文章を表示します (図 1 参照)。

コメントをいくつか

実際、t 検定を必要とする統計問題の種類は複数あります。ここで取り上げた種類の問題は、データ値と各標本データセットとの間に概念的なつながりがないため、対応のない t 検定と呼ばれることがあります。対応のある t 検定と呼ばれる種類の t 検定もあります。対応のある t 検定は、なんらかの指導を行う前のテストの点数と、指導後のテストの点数など、データに前後関係がある場合などに使用します。この場合、点数の各ペアに概念的関係があることになります。

今回説明したウェルチの t 検定は、多くの場合に、一般的なスチューデントの t 検定よりも優れています。スチューデントの t 検定では、一般に 2 つの標本データセットのデータポイント数が同じでなければならず、2 つの標本の分散もほぼ同じであることが求められます。ウェルチの t 検定は、標本サイズが異なっていても機能し、標本の分散が異なる場合も安定しています。

今回説明した種類の t 検定は、両側検定と呼ばれます。これは、程度の差はあるものの、2 群の平均が等しいかどうかを判断することを目標にする問題と同じです。片側 t 検定は、1 つの群の平均がもう 1 つの群の平均よりも大きくなるかどうかを判断することを目標にする状況に使用します。片側 t 検定を実行する場合は、両側 p-value を 2 で除算します。

t 検定の結果の解釈には、慎重さが求められます。t 検定で p-value=0.008 という計算結果が求められた場合、これに基づいて、「男子学生と女子学生の真の母集団の平均点が等しくなる可能性は低い」と結論付けることには妥当性がありますが、「男子学生の平均点は女子学生の平均点よりも高くなる」と結論付ける妥当性は前者よりもはるかに低くなります。

t 検定の代わりとして、マンホイットニーの U 検定と呼ばれる検定もあります。どちらの手法も、標本に基づいて 2 つの母集団の平均が等しいかどうかを推測しますが、マンホイットニーの U 検定では、統計的推測はほとんど行われず、より慎重な結論になります (調査対象の平均が異なるという結論になることはほとんどありません)。

t 検定は、群が 2 つの状況に限定されます。3 つ以上の群の平均を調査する問題では、f 検定という分析手法を使用します。

Dr. James McCaffreyは、ワシントン州レドモンドにある Microsoft Research に勤務しています。これまでに、Internet Explorer、Bing などの複数のマイクロソフト製品にも携わってきました。McCaffrey 博士の連絡先は、jammc@microsoft.com (英語のみ) です。

この記事のレビューに協力してくれた技術スタッフの Kirk Olynyk (Microsoft Research) に心より感謝いたします。