[アーティクル]
07/30/2015

2015 年 7 月

Volume 30 Number 7

テストの実行 - C# を使用した線形回帰

James McCaffrey 線形回帰問題の目標は、1 つ以上の数値予測変数に基づいて数値変数の値を予測することです。たとえば、学歴、就労年数、性別 (男性 = 0、女性 = 1) に基づいて年収を予測するといった考え方です。

予測する変数を通常、従属変数と呼びます。予測の基になる変数を通常、独立変数と呼びます。予測の基になる変数が 1 つしかない場合、この手法を単線形回帰と呼びます。予測の基になる変数が 2 つ以上ある場合、この手法を多重線形回帰 (多変量線形回帰) と呼びます。

今回の主旨を理解するには、図 1 のデモプログラムを見るのが一番です。この C# デモプログラムは、学歴、就労年数、性別に基づいて年収を予測します。デモでは、まず、10 個の合成データ項目を生成します。学歴は 12 ～ 16 の値、就労年数は 10 ～ 30 の値、性別はインジケーター変数で、男性は参照値 0、女性は 1 としてコード化されています。年収は最終列に千ドル単位で示しています。実際のシナリオでは、おそらく、MatrixLoad といった名前のメソッドを使用して、テキストファイルからデータを読み取ることになります。

図 1 C# を使用した線形回帰

合成データを生成後、デモプログラムはこのデータを使用して、いわゆる計画行列を作成します。計画行列は、先頭列すべてに 1.0 という値を追加したデータ行列にすぎません。線形回帰にはさまざまなアルゴリズムを利用できます。手を加えていないデータの行列に使用できるアルゴリズムもあれば、計画行列を使用するアルゴリズムもあります。今回のデモでは、計画行列を必要とする手法を使用します。

計画行列を作成後、デモプログラムは、4 つの係数値 (12.0157, 1.0180, 0.5489, -2.9566) を求めます。これらの係数を b 値またはベータ値と呼ぶこともあります。最初の値 12.0157 を通常インターセプトと呼びます。これは、他の予測変数に関連付けられない定数です。2 番目、3 番目、および 4 番目の係数値 (1.0180, 0.5489, -2.9566) はそれぞれ、学歴、就労年数、性別に関連付けられます。

図 1 の末尾近くの出力は、この係数値を使用して、学歴が 14、就労年数が 12 年、性別が 0 (男) の仮想人物の年収を予測した結果です。予測年収は 32.86 で、以下のように計算しています。

income = 12.0157 + (1.0180)(14) + (0.5489)(12) + (-2.9566)(0)
                = 12.0157 + 14.2520 + 6.5868 + 0
                = 32.86

つまり、線形回帰を使用して予測を行うには、予測値に対応する係数値を乗算し、それらの積を合計します。実にシンプルです。先頭のインターセプト値 (この例では 12.0157) は、常に値 1 を持つ予測変数に対応付けられる係数と考えてもかまいません。このことから、計画行列の先頭にある 1.0 値の列の意味がある程度わかります。

線形回帰問題の本質は、データそのもの、または計画行列を使用して係数の値を計算することにあります。これはそれほど簡単なことではありません。デモでは、最小二乗法とも呼ばれる、閉形式逆行列手法を使用しています。係数値を求める手法には他にも、反復再重み付け最小二乗法、最尤推定法、リッジ回帰、勾配降下などがあります。

図 1 では、予測を行う前に、デモプログラムは R-squared (決定係数) 値というメトリックを計算しています。R-squared は、手を加えていないデータに予測モデルが適合する度合いを示す 0 ～ 1 の値です。これは、「モデルによって示される変化量の割合」と表現されることもあります。大まかに言えば、R-squared が 1 に近づくほど優れた予測モデルになります。デモ値 0.7207 (72%) は、実データに対して比較的高い (良い) 結果と考えられます。

今回は、中級以上のプログラミングスキルがあることを前提にしていますが、線形回帰について何も知らなくても問題ありません。デモプログラムは長すぎてここにすべて掲載することはできませんが、完全なソースコードは本稿付属のコードダウンロードから入手できます。

線形回帰について

線形回帰については、図を見ていただくのが最もわかりやすいでしょう。図 2 のグラフをご覧ください。グラフのデータは、1 つの変数 (就労年数) だけで年収を予測するようすを示しています。それぞれの赤い点は、データポイントに対応します。たとえば、左端のデータ項目は、就労年数 = 10 に対して年収 = 32.06 という値を持ちます。線形回帰では、1 つはインターセプト、もう 1 つは就労変数に対する 2 つの係数を求めます。それぞれの係数は 27.00 と 0.43 であることがわかります。

図 2 独立変数を 1 つ使用する線形回帰

係数値によって、青い直線の方程式が決まります (図 2 参照)。この直線 (係数) は、実際のデータポイント (yi) と予測データポイント (fi) 間の偏差の平方和を最小化します。10 個の偏差のうち 2 個の偏差を破線で示しています (図 2 参照)。図の最初の偏差は yi - fi = 28.6 - 32.6 = -4.0 です。偏差は正の値にも負の値にもなります。偏差を 2 乗しなかった場合、負の値と正の値は互いに打ち消し合います。

図 2 のグラフは、1 つの独立変数だけでシンプルな線形回帰が機能するようすを示しています。同じ考え方を拡張したものが多変量線形回帰です。多変量線形回帰では複数の独立変数を使用して、偏差の平方和を最小化する係数を求めます。

わかりやすく言えば、線形回帰とは、一連のデータポイントに対して最適な直線を見つけることです。この最適な直線を予測に使用します。たとえば、仮想の人物の就労経験が 25 年であれば、青い直線から年収は 38 であると予測されます (図 2 参照)。

最小二乗方程式を解く

線形回帰の問題に n 個の予測変数がある場合、それぞれの予測変数に 1 つずつと、さらにインターセプト値を加えて、n+1 個の係数値を求めなければなりません。デモプログラムでは、最も基本的な手法を使用して係数値を求めています。係数値は、多くの場合、やや難解な方程式を使用して求めます (図 3 参照)。ただ、方程式は見かけほど複雑ではありません。

図 3 行列を使用した線形回帰の係数の解法

ギリシャ文字 β はローマ字の B に似ているため、B を使って係数値を表します。方程式のすべて文字を太字で表しています。文字を太字で表している場合、数学上、その文字は単純なスカラー値 (プレーンな数値) ではなく、複数の値を持つオブジェクト (行列や配列/ベクトル) であることを意味します。大文字の X は計画行列を表します。指数値 T の付いた大文字の X は、計画行列の転置を意味します。* 記号は、行列の乗算です。指数値 -1 は逆行列の意味です。大文字の Y は、従属変数値の列ベクトル (1 列の行列) です。そのため、係数値を求めるには、行列演算を理解しておく必要があります。

図 4 は、行列の転置、行列の乗算、および逆行列を示しています。行列の転置とは、単に行と列を入れ替えることです。たとえば、2 行 3 列の行列があるとします。行列の転置では、元の行列の行が転置行列の列となり、3 行 2 列の行列になります。

図 4 線形回帰の係数を求めるために使用する 3 つの行列演算

行列の乗算を初めて目にすると、わかりにくいと感じるかもしれません。n 行 m 列の行列と m 行 p 列の行列の乗算は、結果が n 行 p 列の行列になります。たとえば、(3 行 4 列の行列) * (4 行 2 列の行列) は 3 行 2 列の行列です。行列の乗算については今回詳しく取り上げませんが、いくつか例を見れば、簡単にプロセスを理解してそのコードを実装できるようになります。

線形回帰の係数値を求めるために必要な 3 つ目の行列演算は、逆行列の計算です。この計算は、理解するのが難しく、実装も困難です。今回の内容を理解するだけであれば、行数と列数が同じ行列 (正方行列) の場合のみ逆行列が定義されることだけを知っておけば十分です。図 4 に示しているのは、3 行 3 列の行列とその逆行列です。

逆行列を求めるアルゴリズムはいくつかあります。デモプログラムでは、ドゥーリトル分解と呼ばれる手法を使用しています。

ここまでの説明をまとめると、予測変数が n 個ある線形回帰問題では n+1 個の係数を求める必要があり、係数を求めるには、行列の転置、行列の乗算、および逆行列を使用する必要があります。転置と乗算は簡単ですが、逆行列を求めるのは困難です。

デモプログラムの構造

デモプログラムを作成するには、Visual Studio を起動して、コンソールアプリケーションプロジェクトテンプレートを選択します。プロジェクトには「LinearRegression」という名前を付けます。このプログラムは .NET Framework との大きな依存関係がないため、どのバージョンの Visual Studio でも機能します。

テンプレートコードがエディターに読み込まれたら、ソリューションエクスプローラーウィンドウで Program.cs ファイルを右クリックし、名前を「LinearRegressionProgram.cs」に変更します。Program クラスの名前が Visual Studio によって自動的に変更されます。エディターウィンドウの上部にある using ステートメントを、最上位レベルの System 名前空間を指定するステートメントを除いてすべて削除します。

デモプログラムの全体構造を図 5 に示します (スペースを節約するために少し編集しています)。プログラムの制御ロジックはすべて Main メソッドに含めています。このデモプログラムでは、オブジェクト指向プログラミング (OOP) のアプローチではなく、静的メソッドのアプローチを使用します。

図 5 線形回帰デモプログラム構造

using System;
namespace LinearRegression
{
  class LinearRegressionProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("Begin linear regression demo");
      // Generate synthetic data
      // Create design matrix
      // Solve for LR coefficients
      // Calculate R-squared value
      // Do a prediction
      Console.WriteLine("End linear regression demo");
      Console.ReadLine();
    }
    static double Income(double x1, double x2,
      double x3, double[] coef) { . . }
    static double RSquared(double[][] data,
      double[] coef) { . . }
    static double[][] DummyData(int rows,
      int seed) { . . }
    static double[][] Design(double[][] data) { . . }
    static double[] Solve(double[][] design) { . . }
    static void ShowMatrix(double[][] m, int dec) { . . }
    static void ShowVector(double[] v, int dec) { . . }
    // ----------
    static double[][] MatrixTranspose(double[][] matrix)
      { . . }
    static double[][] MatrixProduct(double[][] matrixA,
      double[][] matrixB) { . . }
    static double[][] MatrixInverse(double[][] matrix)
      { . . }
    // Other matrix routines here
  }
} // ns

メソッド Income は係数値の配列を使用して、学歴、就労年数、性別の入力パラメーターから予測した年収を返します。メソッド RSquared は、データと係数から、モデルの R-squared 値を返します。メソッド DummyData は、デモで使用する合成データを生成します。

メソッド Design はデータの行列を受け取り、先頭列に 1.0 値を持つように拡張した計画行列を返します。メソッド Solve は計画行列を受け取り、行列演算を使用して線形回帰の係数を求めます。

難しい作業の大部分は、行列演算を実行する一連の静的メソッドで処理します。デモプログラムは、可能な限り最もシンプルな方法 (配列の配列) で行列を定義しています。他にもプログラム定義の Matrix クラスを作成する方法もありますが、個人的な意見としては無駄に複雑になると感じています。場合によっては、通常の配列の方がプログラム定義のオブジェクトよりも優れています。

メソッド MatrixTranspose は、転置行列を返します。メソッド MatrixProduct は、2 つの行列を乗算した結果を返します。メソッド MatrixInverse は、逆行列を返します。このデモには多くのヘルパーメソッドを用意しています。特に、メソッド MatrixInverse では、ヘルパーメソッドとして MatrixDuplicate、MatrixDecompose、および HelperSolve を呼び出します。

Solve メソッド

線形回帰デモプログラムの心臓部はメソッド Solve です。このメソッドの定義は、次のように始めます。

static double[] Solve(double[][] design)
{
  int rows = design.Length;
  int cols = data[0].Length;
  double[][] X = MatrixCreate(rows, cols - 1);
  double[][] Y = MatrixCreate(rows, 1);
...

入力パラメーターは 1 つだけで、計画行列を受け取ります。他にも、ソースデータ行列を渡してから、Solve メソッドでヘルパーメソッド Design を呼び出して計画行列を取得するようにする方法もあります。ヘルパーメソッド MatrixCreate は、指定された数の行と列を持つ行列用の領域の確保し、その行列を返します。ローカルの行列 X は、独立変数の値を保持します (先頭には 1.0 値を保持します)。ローカルの行列 Y には 1 列しかなく、従属変数 (このデモでは年収) の値を保持します。

次に、行列 X と行列 Y の各セルに、計画行列の値を使用して値を設定します。

int j;
for (int i = 0; i < rows; ++i)
{
  for (j = 0; j < cols - 1; ++j)
  {
    X[i][j] = design[i][j];
  }
  Y[i][0] = design[i][j]; // Last column
}

インデックス変数 j は、入れ子になった for ループの外で宣言しているため、行列 Y への値の設定に使用できます。行列 X と行列 Y を使用して、線形回帰の係数を求めることができます (図 3 の方程式を参照)。

...
  double[][] Xt = MatrixTranspose(X);
  double[][] XtX = MatrixProduct(Xt, X);
  double[][] inv = MatrixInverse(XtX);
  double[][] invXt = MatrixProduct(inv, Xt);
  double[][] mResult = MatrixProduct(invXt, Y);
  double[] result = MatrixToVector(mResult);
  return result;
} // Solve

このデモでは、行列 X のサイズは 10 行 4 列なので、転置 Xt は 4 行 10 列になります。Xt と X の積は 4 行 4 列になり、逆行列 inv も 4 行 4 列になります。一般的に、n 個の独立した予測変数を持つ線形回帰問題で逆行列手法を使用する場合、(n+1) 行 (n+1) 列の逆行列を求める必要があります。つまり、逆行列の手法は、予測変数の数が多い線形回帰問題には適していません。

4 行 4 列の逆行列と 4 行 10 列の転置行列 (コードの invXt) の積は 4 行10 列の行列です。コードの invXt と 10 行 1 列の行列 Y との積 mResult は、4 行 1 列の行列になりますす。これらの値が求める係数です。利便性を目的として、単一列の行列 Y の値は、ヘルパーメソッド MatrixToVector を使用して通常の配列に転送します。

R-Squared の計算

前述のように、メトリック R-squared は、実際のデータポイントと計算した回帰直線の適合度を示す尺度です。数学的に表現すると、R-squared は R2 = 1 - (SSres / SStot) と定義されます。項 SSres は、通常「残差平方和」と呼びます。これは、実際の Y 値と予測した Y 値の差を 2 乗して合計したものです (図 2 参照)。項目 SStot は「総平方和」です。これは、実際の各 Y 値と、すべての Y 値の平均との差を 2 乗して合計したものです。

線形回帰のメトリック R-squared は決定係数とも呼ばれ、r-squared (小文字) という静的メトリックと関係はありますが、別物です。R-squared の解釈はやや複雑で、研究途上の特定の問題領域に依存します。扱いにくいデータや不完全なデータが収集される傾向がある自然科学や社会科学の分野では、0.6 以上の R-squared 値はかなり優れていると考えられます。

調整済み R-squared と呼ばれる回帰モデルによって表される変位の尺度もあります。このメトリックは、予測変数の数とデータ項目の数を考慮に入れます。線形回帰モデルの予測品質という点では、ほとんどの場面でプレーンな R-squared 値で十分です。

まとめ

プログラミング言語を使用して線形回帰を実行する方法の例をインターネットで検索しても、あまり資料は見つかりません。情報が比較的少ない主な理由は 2 つあると考えられます。1 つは、行列演算を使用して線形回帰の係数を求めるのがかなり難しいことです (逆行列の演算が主な原因)。ある意味では、このデモプログラムの MatrixInverse メソッドが、これまで作成した中で最も複雑なコードルーチンかもしれません。もう 1 つは、Excel のデータ分析アドインを追加した Excel ワークシートのプログラムなど、線形回帰を実行できるスタンドアロンのツールが既に豊富に存在するためです。線形回帰の解法をソフトウェアシステムに直接組み込む必要が生じることはあまりありません。

線形回帰は 10 年以上研究されているため、線形回帰の手法を拡張する方法はたくさんあります。たとえば、2 つ以上の予測変数を組み合わせる、相互作用効果という手法を導入することができます。これらの拡張手法は、基本形式の線形回帰と区別するため、一般線形モデルと呼ばれることがあります。

線形回帰は、古典統計学の手法における「Hello World」のような存在だと考えられます。古典統計学と機械学習の間には、明確かつ普遍的で合意済みの区別はありませんが、個人的には、1900 年代初頭に初めて数学者が研究し始めた統計学を古典統計学と考えます。ニューラルネットワーク分類のような機械学習手法は、1950 年代に初めて発表された比較的新しい手法だという考えです。古典統計学の線形回帰は、「テストの実行」コラムのトピックとして何回も取り上げたことのあるロジスティック回帰という機械学習手法に密接に関連しています。

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

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