7. 単回帰分析と重回帰分析

本章では、基礎的な機械学習手法として代表的な単回帰分析重回帰分析の仕組みを、数式を用いて説明します。 また次章では、本章で紹介した数式を Python によるプログラミングで実装する例も紹介します。本章と次章を通じて、数学とプログラミングの結びつきを体験して理解することができます。

本チュートリアルの主題であるディープラーニングの前に、単回帰分析と重回帰分析を紹介することには 2 つの理由があります。 1 つ目は、単回帰分析と重回帰分析の数学がニューラルネットワーク含めたディープラーニングの数学の基礎となるためです。 2 つ目は、単回帰分析のアルゴリズムを通して微分、重回帰分析のアルゴリズムを通して線形代数に関する理解を深めることができるためです。

機械学習手法は、教師あり学習 (supervised learning)教師なし学習 (unsupervised learning)強化学習 (reinforcement learning) に大別され、単回帰分析は教師あり学習に含まれます。 本チュートリアルで扱う多くの手法が教師あり学習です。

教師あり学習の中でも典型的な問題設定は 2 つに大別されます。 与えられた入力変数から、100.1 といった実数値を予測する 回帰 (regression) と、「赤ワイン」、「白ワイン」といったカテゴリを予測する分類 (classification) の 2 つです。

単回帰分析は回帰を行うための手法であり、1 つの入力変数から 1 つの出力変数を予測します。 それに対し、重回帰分析は、複数の入力変数から 1 つの出力変数を予測します。 この両手法は教師あり学習であるため、訓練の際には、入力変数 x と目的変数 t がペアで準備されている必要があります。

回帰分析を行うアルゴリズムでは、以下の 3 ステップを順番に考えていきます。

  • Step 1 : モデルを決める

  • Step 2 : 目的関数を決める

  • Step 3 : 最適なパラメータを求める

7.1. 単回帰分析

7.1.1. 問題設定(単回帰分析)

単回帰分析では、 1 つの入力変数から 1 つの出力変数を予測します。 今回は身近な例として、部屋の広さ x から家賃 y を予測する問題を考えてみます。

7.1.2. Step 1:モデルを決める(単回帰分析)

まずはじめに、入力変数 x と出力変数 y との関係をどのように定式化するかを決定します。 この定式化したものを モデル もしくは 数理モデル と呼びます。

単回帰分析におけるモデルを具体的に考えていきましょう。 例えば、家賃と部屋の広さの組で表されるデータを 3 つ集め、「家賃」を y 軸に、「部屋の広さ」を x 軸にとってそれらをプロットしたとき、次のようになっていたとします。

部屋の広さと家賃

この場合、部屋が広くなるほど、家賃が高くなるという関係が予想されます。 また、この 2 変数間の関係性は直線によって表現を行うことができそうだと考えられます。 そこで、2 つのパラメータ wb によって特徴づけられる直線の方程式

y=wx+b

によって、部屋の広さと家賃の関係を表すことを考えます。 ここで、w重み (weight)bバイアス (bias) の頭文字を採用しています。

部屋の広さと家賃の関係

単回帰分析では、このようにモデルとして直線 y=wx+b を用います。 そして、2 つのパラメータ wb を、直線がデータによくフィットするように調整します。

パラメータで特徴づけられたモデルを用いる場合、与えられた データセット に適合するように最適なパラメータを求めることが目標となります。 今回はデータセットとして部屋の広さ x と家賃 t の組からなるデータの集合を用います。 全部で N 個のデータがあり、n 番目のデータが (xn,tn) と表されるとき、データセットは

D={(x1,t1),(x2,t2),,(xN,tN)}={(xn,tn)}Nn=1

と表すことができます(注釈1)。 これを用いて、新しい x を入力すると、それに対応する t を予測するモデルを訓練します。

7.1.3. 前処理

次のステップに進む前に、データの前処理 (preprocessing) をひとつ紹介します。 データの 中心化 (centering) は、データの平均値が 0 になるように全てのデータを平行移動する処理を言います。 下図は、データ集合 (xn,yn) (n=1,,11) を、平均が (0,0) になるように平行移動する例です。

データの中心化

中心化によるデータ前処理の利点の一つとして、調整すべきパラメータを削減できることが挙げられます。中心化処理を行うことで切片を考慮する必要がなくなるため、データ間の関係性を表現する直線の方程式を、 yc=wxc のように、簡潔に表現可能となります。

中心化によるパラメータ削減

データセット内の入力変数と目的変数の平均をそれぞれ ˉx, ˉt としたとき、中心化後の入力変数と目的変数は、

xc=xˉxtc=tˉt

となります。

以降は記述を簡単にするため、c という添え字を省略し、事前に中心化を行っている前提でデータを扱います。 また、そのデータにフィットさせたいモデルは、

y=wx

と、こちらも添え字 c を省略して説明を行います。

7.1.4. Step 2:目的関数を決める(単回帰分析)

ここでの目標は、部屋の広さと家賃の関係を直線の方程式によってモデル化することです。 このために、予め収集されたいくつかのデータセットを使って、モデルが部屋の広さの値から予測した家賃(予測値)と、その部屋の広さに対応する実際の家賃(目標値)の差が小さくなるように、モデルのパラメータを決定します。

今回は、目的関数として こちらの章 ですでに紹介した予測値と目標値との二乗和誤差 (sum-of-squares error) を用います。 二乗和誤差が 0 のとき、またその時のみ予測値 y は目標値 t と完全に一致(t=y)しています。 n 個目のデータの部屋の広さ xn が与えられたときのモデルの予測値 yn と、対応する目標値 tn との間の二乗誤差は、

(tnyn)2

となります。これを全データに渡って合計したものが以下の二乗和誤差です。

L=(t1y1)2+(t2y2)2++(tNyN)2=Nn=1(tnyn)2

今回用いるモデルは

yn=wxn

であるため、目的関数は

L=Nn=1(tnwxn)2

と書くこともできます。

7.1.5. Step 3:最適なパラメータを求める(単回帰分析)

この目的関数を最小化するようなパラメータを求めます。 ここで、目的関数は差の二乗和であり、常に正の値または 0 をとるような、下に凸な二次関数となっています。 (一般的には多くの場合において、最適なパラメータを用いてもモデルがすべてのデータを完全に表現できず、目的関数の値は 0 にはなりません。)

目的関数の値が最小となる点を求める際には、微分の知識が有用です。 微分では、対象とする関数の接線の傾きを求めることができます。凸関数では、この接線の傾きが 0 である点において、関数の最小値、もしくは最大値が得られます。 今回は、目的関数が x に関する二次関数となっているため、下図のように重み w に関する接線の傾きが 0 であるときに、目的関数の値が最小となります。

目的関数を最小化するパラメータ

それでは、具体的に今回定めた目的関数 L をパラメータである w で微分してみましょう。 微分に関する基本的な計算や性質は こちらの章で紹介しました。

wL=wNn=1(tnwxn)2

ここで、微分の線形性から、和の微分は、微分の和となるため、

wL=Nn=1w(tnwxn)2

と変形できます。

次に、総和()の中の各項に着目すると、

w(tnwxn)2

となっており、この部分は tnwxn()2 という関数の合成関数になっています。 そこで、un=tnwxnf(un)=u2n とおいて計算すると、

w(tnwxn)2=wf(un)=unwf(un)un=xn(2un)=2xn(tnwxn)

が得られます。

これを L/w の式に戻すと、

wL=Nn=1w(tnwxn)2=Nn=12xn(tnwxn)

となります。

この導関数の値が 0 となるときの w が、目的関数を最小にするパラメータです。 そこで、wL=0 とおいてこれを w について解きます。

wL=02Nn=1xn(tnwxn)=02Nn=1xntn+2Nn=1wx2n=02Nn=1xntn+2wNn=1x2n=0wNn=1x2n=Nn=1xntn

以上より、

w=Nn=1xntnNn=1x2n

と求まりました。これを最適な w と呼びます。 この値は、与えられたデータセット D={xn,tn}Nn=1 のみから決定されています。

7.1.6. 数値例

例題にあげていた数値でパラメータ w を求めてみましょう。

まずはデータの中心化を行うために、平均の値を事前に算出します。

ˉx=13(20+40+60)=40ˉt=13(60000+115000+155000)=110000

この平均の値を使い、全変数に対して中心化を行うと、

x1=2040=20x2=4040=0x3=6040=20t1=60000110000=50000t2=115000110000=5000t3=155000110000=45000

となります。

これらの中心化後の値を用いて、最適なパラメータ w を計算すると、

w=Nn=1xntnNn=1x2n=x1t1+x2t2+x3t3x21+x22+x23=20×(50000)+0×5000+20×45000(20)2+02+202=2375

と求まります。 したがって、家賃が 1 m2 増えるごとに、2375 円家賃が上昇しているとわかります。

この w を用いて決定される直線 y=2375x と、学習データとして用いた 3 つの点をプロットした図が以下です。

家賃の予測値と目標値(中心化後)

この直線上の点の y の値が、対応する x の値に対するここで訓練したモデルによる予測値です。 ここで、x 軸で負の値をとっていますが、これは中心化後であることに注意してください。

訓練済みのモデルを使って新しいサンプル xq に対する家賃の予測を行います。 例えば、部屋の広さが 50 m2 の場合の家賃に対する予測値を計算する推論を行いましょう。

yc=wxcyqˉt=w(xqˉx)yq=w(xqˉx)+ˉt=2375(5040)+110000=133750

このように、部屋の広さが 50 m2 の場合の家賃が 133,750 円であると予測することができました。 上記のように、モデル化の際に中心化を行っていた処理を推論の際には元に戻して計算しましょう。

7.2. 重回帰分析

7.2.1. 問題設定(重回帰分析)

重回帰分析も単回帰分析の場合と同様に家賃を予測する問題を題材にして説明します。 単回帰分析の場合と異なり、入力変数として「部屋の広さ」だけでなく、「駅からの距離」や「犯罪発生率」などの変数を合わせて考慮することにします。 部屋の広さを x1、駅からの距離を x2、…、犯罪発生率を xM といった形で表し、M 個の入力変数を扱うことを考えてみましょう。

7.2.2. Step 1:モデルを決める(重回帰分析)

単回帰分析では、

y=wx+b

のように直線の方程式をモデルとして用いていました。重回帰分析でも、

y=w1x1+w2x2++wMxM+b

のように、単回帰分析と似た形のモデルを定義します。 単回帰分析の際は二次元平面を考え、その平面上に存在するデータに最もよくフィットする直線を考えましたが、今回は M 次元空間に存在するデータに最もよくフィットする 直線を考えることになります。

重回帰分析のモデルは総和の記号を使って表記すると、

y=Mm=1wmxm+b

と書くことができます。

ここでバイアス b の扱い方を改めて考えます。 単回帰分析では、中心化を前処理として施し、バイアス b を省略することで、簡潔に定式化することができました。

重回帰分析では、 M 個の重み w1,w2,,wM と 1 個のバイアス b があり、合わせて M+1 個のパラメータが存在します。これらのパラメータをうまく定式化することを考えます。 そこで、今回は x0=1w0=b とおくことで、

y=w1x1+w2x2++wMxM+b=w1x1+w2x2++wMxM+w0x0=w0x0+w1x1++wMxM=Mm=0wmxm

のように b を総和の内側の項に含めて、簡潔に表記できるようにします。 (ここで、 記号の下部が m=1 ではなく m=0 となっていることに注意してください。)

さらに、ここから線形代数で学んだ知識を活かして、数式を整理していきます。 上式をベクトルの内積を用いて表記しなおすと、

y=w0x0+w1x1++wMxM=[w0w1wM][x0x1xM]=wTx

のように、シンプルな形式で表現することができます。 このモデルが持つパラメータは前述の通り M+1 個あり、M+1 次元のベクトル w で表されています。 重回帰分析では、この w のすべての要素について最適な値を求めます。

7.2.3. Step 2:目的関数を決める(重回帰分析)

単回帰分析の例と比べると、入力変数は増えましたが、家賃を目標値としている点は変わっていません。 そこで、単回帰分析と同じ目的関数

L=(t1y1)2+(t2y2)2++(tNyN)2

を用います。 この目的関数は、ベクトルの内積を用いて表記し直すと、

L=(t1y1)2+(t2y2)2++(tNyN)2=[t1y1t2y2tNyN][t1y1t2y2tNyN]=(ty)T(ty)

と書くことができます。

ここで、内積には交換法則が成り立つため、wTxxTw と書くこともできます。これを利用して、モデルの方程式 y=wTx を、以下のように変形します。

y=[y1y2yN]=[xT1wxT2wxTNw]=[xT1xT2xTN]w

さらに、xTn=[xn0, xn1, xn2, , xnM] (n=1,,N) と展開すると、

y=[x10x11x12x1Mx20x21x22x2MxN0xN1xN2xNM][w0w1w2wM]=Xw

と表記できます。 ここで、N×M 行列 X は、各行が各データを表しており、各列が各入力変数を表しています。 このような行列を、デザイン行列 (design matrix) と呼びます。 各列はそれぞれ入力変数の種類に対応しており、例えば、部屋の広さや、駅からの距離などです。

各行が表すデータ点がどのように表されているかを説明するため、具体的な数値例を挙げます。 部屋の広さ =50m2 、駅からの距離 =600m 、犯罪発生率 =2% という 3 つの入力変数を考える場合、M=3 であり、これが n 個目のデータのとき、xTn は、

xTn=[1506000.02]

となります。先頭に 1 があるのは、Step 1 で x0=1 と定めたためです。

7.2.4. Step 3:パラメータを最適化する(重回帰分析)

それでは、目的関数 L を最小化するモデルのパラメータベクトル w を求めましょう。 単回帰分析と同様に、目的関数をパラメータで微分して 0 とおき、w について解きます。

まずは目的関数に登場している予測値 y を、パラメータ w を用いた表記に置き換えます。

L=(ty)T(ty)=(tXw)T(tXw)={tT(Xw)T}(tXw)=(tTwTXT)(tXw)

ここで、転置の公式 (AB)T=BTAT を用いました。

さらに分配法則を使って展開すると、

L=tTttTXwwTXTt+wTXTXw

となります。この目的関数に対しパラメータの w についての偏微分を計算します。

その前に、この式をもう少し整理します。 まず、(1)T=1 のように、スカラは転置しても変化しません。 また、wRM+1XRN×(M+1) であり、XwRN となることから、tRN との間の内積 tTXwR は、スカラになります。 したがって、

(tTXw)T=tTXw

が成り立ちます。

さらに、転置の公式 (ABC)T=CTBTAT より、

(tTXw)T=wTXTt

も成り立ちます。以上から、

(tTXw)T=tTXw=wTXTt

が導かれます。目的関数 L は、この式を利用して、

L=tTt2tTXw+wTXTXw

と変形できます。

ここで、w に関する偏微分を行いやすくするため、w 以外の定数項を一つにまとめます。

L=tTt2tTXw+wTXTXw=tTt2(XTt)Tw+wTXTXw=c+bTw+wTAw

すると、w に関する二次形式で表現できました。 ここで、

A=XTXb=2XTtc=tTt

とおいていることに注意してください。

それでは、目的関数を最小にするパラメータ w の求め方を考えます。 目的関数はパラメータ w に関して二次関数になっています。 まずは、w 以外のベクトルや行列に、具体的な数値を当てはめて考えてみましょう。 例えば、

w=[w1w2],A=[1234],b=[12],c=1

とおきます。すると、目的関数の値は

L=wTAw+bTw+c=[w1w2][1234][w1w2]+[12][w1w2]+1=[w1w2][w1+2w23w1+4w2]+w1+2w2+1=w1(w1+2w2)+w2(3w1+4w2)+w1+2w2+1=w21+5w1w2+4w22+w1+2w2+1

となります。これを w1,w2 に関して整理すると、

L=w21+(5w2+1)w1+(4w22+2w2+1)=4w22+(5w1+2)w2+(w21+w1+1)

となり、w1,w2 それぞれについて二次関数になっていることが分かります。

目的関数 Lw1 の二次関数、w2 の二次関数と見たとき、L は、下図のような概形となっています。

目的関数の概形(2 次元)

さらに、各次元が w1,w2,L を表す 3 次元空間上においては、 L の概形は下図のようになっています。

目的関数の概形(3 次元)

ここでは 2 つのパラメータ w1w2 について図示していますが、目的関数が 任意の M 個の変数 w0,w1,w2,,wM によって特徴づけられている場合でも、目的関数がそれぞれのパラメータについて二次形式になっている限り、同様のことが言えます。 すなわち、M+1 個の連立方程式、

{w0L=0w1L=0     wML=0

を解けば良いということになります。 これはベクトルによる微分を用いて表記すると、以下のようになります。

[w0Lw1LwML]=[000]wL=0

上式を w について解くために、以下のような式変形を行います。 式変形の途中で理解できない部分があった場合は、こちらの章 を読み返してみてください。 まずは、左辺について整理を行います。

wL=w(c+bTw+wTAw)=w(c)+w(bTw)+w(wTAw)=0+bT+wT(A+AT)

これを 0 とおき、Ab を展開すると

2(XTt)T+wT{XTX+(XTX)T}=02tTX+2wTXTX=0wTXTX=tTX

のように式変形できます。

ここで両辺を転置して、 wT の転置を直しておきましょう。

(wTXTX)T=(tTX)TXTXw=XTt

ここで、XTX逆行列が存在すると仮定して、両辺に左側から (XTX)1 を掛けると、

(XTX)1XTXw=(XTX)1XTtIw=(XTX)1XTtw=(XTX)1XTt

が導かれます。特に、この最後の式を正規方程式 (normal equation) と呼びます。 上式は、与えられたデータを並べたデザイン行列 X と、各データの目標値を並べたベクトル t から、最適なパラメータ w を計算しています。 I は単位行列を表します。

w を求める際に気をつけたいこととして、以下の誤った式変形をしてしまう例が挙げられます。

XTXw=XTt(XT)1XTXw=(XT)1XTtXw=tX1Xw=X1tw=X1t

上記の式変形は一般には成立しません。 この式変形が可能かどうかは、X1 が存在するかどうか、に関わっています。 サンプルサイズ N と独立変数の数 M+1 が等しくない場合、XRN×(M+1)正方行列ではないため、逆行列 X1 を持ちません。 従って、上式の 2 行目の変形を行うことはできません(逆行列が求まるためのより厳密な条件についてはここでは省略します)。

一方、 XTX(M+1)×(M+1) 行列であり、その形はサンプルサイズ N に依存することなく、常に正方行列となるため、これを利用して式変形を行います。

新しい入力変数の値 xq=[x1,,xM]T に対して、対応する目標値 yq を予測するためには、訓練により決定されたパラメータ w を用いて、

yq=wTxq

のように計算します。


注釈 1

データセット中の各データのことをデータ点(datum)ということがあります。具体的には、上の説明で同上した D 中の各 (x1,t1) などのことです。

▲上へ戻る