マンションの家賃は何で決まるだろうか。部屋が広いほど家賃は高く、駅から遠いほど安くなりそうだ。このように、1つの量（家賃）を、複数の量（広さや駅からの距離など）を使って説明したり予測したりする手法を重回帰分析（multiple regression analysis）という。説明に使う量が1つだけの場合は単回帰分析と呼ばれ、データに直線をあてはめる。重回帰分析はこれを複数の量に広げたものであり、売上の予測、病気のリスク要因の分析、経済データの分析など、データ分析の最も基本的な道具として広く使われている。
重回帰分析の土台は最小二乗法（method of least squares）である。最小二乗法を初めて出版したのはフランスの数学者アドリアン＝マリ・ルジャンドル（Adrien-Marie Legendre）で、1805年に出した彗星の軌道を決める方法についての著作の付録で、この方法を明快に述べた。一方、ドイツの数学者カール・フリードリヒ・ガウス（Carl Friedrich Gauss）は1809年に天体の運動についての著作で最小二乗法を発表し、1795年からこの方法を使っていたと主張したため、2人の間で先取権の争いが起こった。ガウスの名を一躍有名にしたのが、ケレス（現在は準惑星に分類される天体）の軌道の計算である。1801年1月1日にイタリアの天文学者ジュゼッペ・ピアッツィ（Giuseppe Piazzi）が発見したケレスは、40日ほど追跡されたのち太陽の光にまぎれて見えなくなった。当時24歳のガウスは最小二乗法を使って限られた観測から軌道を計算し、その予測をもとに、ケレスは1801年12月に再び見つけ出された。
「回帰」という名前は、イギリスの科学者フランシス・ゴルトン（Francis Galton）に由来する。ゴルトンは1886年の論文「遺伝における身長の平凡への回帰」（Regression towards Mediocrity in Hereditary Stature）で、205組の両親とその成人した子どもたちの身長を調べた（女性の身長は1.08倍して男性と比べられるようにした）。その結果、背の高い両親の子どもは平均より背が高いものの両親ほどではなく、子どもの平均からのずれは両親のずれのおよそ3分の2に縮むことを見いだした。平均へ戻っていくこの現象を、ゴルトンは regression（回帰）と呼んだ。もともとは生物の遺伝の現象を指す言葉だったが、その後ユール（Udny Yule）やピアソン（Karl Pearson）によってより一般的な統計の枠組みへと広げられ、現在の回帰分析につながった。2つの変数が正規分布に従う場合の回帰直線については、多変量正規分布の記事でも扱っている。
本記事では、説明に使う量が1つの単回帰から始めて重回帰へと広げ、最小二乗推定量の求め方（正規方程式）、偏回帰係数の意味、行列を使った表し方と射影としての見方、決定係数と自由度調整済み決定係数、最小二乗推定量の性質（不偏性・分散・最小分散性）、多項式回帰までを解説する。回帰係数の検定は別の記事で扱う。
単回帰から重回帰へ
単回帰：直線をあてはめる
具体例から始めよう。ある駅の周辺にある6つの賃貸物件について、広さ、駅からの徒歩時間、家賃を調べたところ、次のようになったとする。
| 物件 | 広さ x_1（㎡） | 駅徒歩 x_2（分） | 家賃 y（万円） |
|---|---|---|---|
| A | 22 | 7 | 7.0 |
| B | 26 | 11 | 6.8 |
| C | 28 | 8 | 8.7 |
| D | 32 | 12 | 8.1 |
| E | 34 | 7 | 9.5 |
| F | 38 | 15 | 7.9 |
まず広さだけで家賃を説明してみる。広さを x、家賃を y として
という式を考える。\beta_0 + \beta_1 x は直線で、\varepsilon はその直線からのずれ（誤差）を表す。このモデルを単回帰モデルという。切片 \beta_0 と傾き \beta_1 をデータから決めるには、図1のように各点と直線の縦のずれ（残差）を考え、残差の2乗の合計が最も小さくなる直線を選ぶ。これが最小二乗法である。
このデータであてはめた直線は \hat{y} = 5.0 + 0.1x となる（求め方は後で説明する）。\hat{y}（ワイハット）は「直線から予測した家賃」を表す。傾きの0.1は「広さが1㎡大きいと、家賃が平均0.1万円（1000円）高い」ことを意味するが、点は直線のまわりに大きくばらついている。たとえば E と F はどちらも広いが、F は家賃が安い。F は駅から徒歩15分と遠いからだろう。
説明変数を増やす
そこで駅からの徒歩時間も使い、広さを x_1、駅徒歩を x_2 として
と考える。x_1 と x_2 を2本の横軸、y を高さにとると、\beta_0 + \beta_1 x_1 + \beta_2 x_2 は3次元空間の中の平面を表す。単回帰が点の集まりに直線をあてはめたのに対し、説明に使う量が2つの重回帰では平面をあてはめることになる。説明に使う量がさらに増えると図には描けないが、考え方は同じである。
一般に、予測したい量 y を目的変数（被説明変数）、説明に使う量 x_1, \dots, x_d を説明変数という。
n 個のデータ (x_{i1}, \dots, x_{id}, y_i)（i = 1, \dots, n）について
が成り立つとするモデルを重回帰モデルという。\varepsilon_1, \dots, \varepsilon_n は互いに独立に正規分布 N(0, \sigma^2) に従う誤差である。
\beta_1, \dots, \beta_d を回帰係数（重回帰では特に偏回帰係数）、\beta_0 を切片という。誤差の分散 \sigma^2 も未知のパラメータである。誤差が正規分布に従うという仮定は、後で述べる性質のうち「最小分散性」の一部や、回帰係数の検定で使う。最小二乗推定量の求め方や不偏性などは、誤差の平均が0、分散が \sigma^2 で、互いに無相関でありさえすれば成り立つ。
最小二乗推定量
残差平方和を最小にする
説明変数が2つの場合で考える。係数をある値 \beta_0, \beta_1, \beta_2 に決めたとき、i 番目のデータの予測値は \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2}、実際の値とのずれは y_i - (\beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2}) である。これらの2乗の合計
を残差平方和という。S を最小にする係数を最小二乗推定量（least squares estimator）といい、\hat{\beta}_0, \hat{\beta}_1, \hat{\beta}_2 と書く。
S は \beta_0, \beta_1, \beta_2 の2次式なので、最小になる点では、どの変数で偏微分しても0になる。たとえば \beta_1 で偏微分すると、合成関数の微分により
となる。\beta_0、\beta_2 についても同様に計算する（\beta_0 では -x_{i1} の代わりに -1 が掛かる）。最小二乗推定量での残差を e_i = y_i - \hat{\beta}_0 - \hat{\beta}_1 x_{i1} - \hat{\beta}_2 x_{i2} と書くと、3つの偏微分が0になる条件は次のように表せる。
これらを正規方程式（normal equation）という。「残差の合計が0」「残差と各説明変数の積の合計が0」という、覚えやすい形をしている。
平均からの偏差で書き直す
1つ目の式 \sum e_i = 0 を書き下して両辺を n で割ると
が得られる（\bar{x}_1, \bar{x}_2, \bar{y} はそれぞれの平均）。これはあてはめた平面が平均の点 (\bar{x}_1, \bar{x}_2, \bar{y}) を通ることを意味する。切片は \hat{\beta}_0 = \bar{y} - \hat{\beta}_1 \bar{x}_1 - \hat{\beta}_2 \bar{x}_2 と表せるので、これを残差の式に代入すると
となり、残差は平均からの偏差だけで書ける。ここで、偏差の2乗和や積和を次のようにおく。
次に、2つ目の式 \sum x_{i1} e_i = 0 も偏差だけの形に直したい。そのために、x_{i1} を偏差 x_{i1} - \bar{x}_1 に取りかえても、残差との積の和は変わらないことを確かめておく。
1行目は、かっこを外して2つの和に分けただけである。2行目では、平均 \bar{x}_1 が i によらない同じ数なので、\bar{x}_1 e_1 + \bar{x}_1 e_2 + \cdots + \bar{x}_1 e_n = \bar{x}_1 (e_1 + e_2 + \cdots + e_n) のように和の外にくくり出した。3行目では、1つ目の正規方程式 \sum e_i = 0 を使った。e_i は最小二乗推定量での残差であり、最小二乗推定量は3つの正規方程式をすべて同時に満たすので、このように1つ目の式を使って2つ目の式を書きかえてよい。連立方程式を解くときに、ある式からほかの式の何倍かを引くのと同じ操作で、解は変わらない（ここでは、2つ目の式から1つ目の式の \bar{x}_1 倍を引いたことになる）。つまり、残差の合計が0なので、x_{i1} から同じ数を引いても、残差との積の和は変わらないのである。
合計が0になる3つの数 (e_1, e_2, e_3) = (1, -3, 2) と、(x_1, x_2, x_3) = (4, 5, 9)（平均6）で比べると
となり、どちらも7である。2つの差は 6 \times (1 - 3 + 2) = 6 \times 0 = 0 になっている。
したがって、2つ目の正規方程式 \sum x_{i1} e_i = 0 は \sum (x_{i1} - \bar{x}_1) e_i = 0 と同じ式である。ここに、偏差で書いた残差の式 e_i = (y_i - \bar{y}) - \hat{\beta}_1 (x_{i1} - \bar{x}_1) - \hat{\beta}_2 (x_{i2} - \bar{x}_2) を代入して展開すると（和はすべて i = 1 から n まで）
となる。最後の行では、3つの和がそれぞれ S_{1y}、S_{11}、S_{12} の定義そのものであることを使った。平均を引いた形に直したことで、切片 \hat{\beta}_0 が式から消え、未知数は \hat{\beta}_1 と \hat{\beta}_2 の2つだけになった。3つ目の式 \sum x_{i2} e_i = 0 も、同じように x_{i2} を x_{i2} - \bar{x}_2 に取りかえて展開すると 0 = S_{2y} - \hat{\beta}_1 S_{12} - \hat{\beta}_2 S_{22} となる。よって、偏回帰係数は次の連立方程式で決まる。
この連立方程式を解くと
が得られる。
単回帰の場合
説明変数が x_1 の1つだけの単回帰 y = \beta_0 + \beta_1 x_1 + \varepsilon でも、手順はまったく同じである。残差は e_i = y_i - \hat{\beta}_0 - \hat{\beta}_1 x_{i1} で、正規方程式は \sum e_i = 0 と \sum x_{i1} e_i = 0 の2つになる。1つ目の式から \hat{\beta}_0 = \bar{y} - \hat{\beta}_1 \bar{x}_1 となるので、残差は e_i = (y_i - \bar{y}) - \hat{\beta}_1 (x_{i1} - \bar{x}_1) と偏差だけで書ける。2つ目の式も、x_{i1} を x_{i1} - \bar{x}_1 に取りかえてから（1つ目の式 \sum e_i = 0 より値は変わらない）この残差を代入すると
となる。重回帰の連立方程式 S_{11} \hat{\beta}_1 + S_{12} \hat{\beta}_2 = S_{1y} から、\hat{\beta}_2 を含む項がなくなった形である。これを解いて、次の公式が得られる。
分子と分母を同じ n - 1 で割ると、分子は x_1 と y の標本共分散、分母は x_1 の標本分散になる。つまり単回帰の傾きは「x_1 と y の共分散 ÷ x_1 の分散」である。共分散は、x_1 が平均より大きいときに y も平均より大きくなりやすいかを表す量で、単位は（x_1 の単位）×（y の単位）である。これを x_1 の分散（単位は x_1 の単位の2乗）で割ると、「x_1 が1増えると y が平均どれだけ増えるか」という傾きの単位になる。2変数の正規分布で X_1 から X_2 を予測する回帰直線の傾き \rho \dfrac{\sigma_2}{\sigma_1} = \dfrac{\sigma_{12}}{\sigma_1^2} も、共分散を分散で割った同じ形をしている。
計算例：家賃のデータ
冒頭の6物件のデータで、実際に偏回帰係数を計算してみよう。データをもう一度示す。
| 物件 | 広さ x_1（㎡） | 駅徒歩 x_2（分） | 家賃 y（万円） |
|---|---|---|---|
| A | 22 | 7 | 7.0 |
| B | 26 | 11 | 6.8 |
| C | 28 | 8 | 8.7 |
| D | 32 | 12 | 8.1 |
| E | 34 | 7 | 9.5 |
| F | 38 | 15 | 7.9 |
手順1（平均を求める）
手順2（偏差の2乗和と積和を求める） 各物件の平均からの偏差を a_i = x_{i1} - \bar{x}_1、b_i = x_{i2} - \bar{x}_2、c_i = y_i - \bar{y} とおき、それらの2乗や積を表にまとめる。
| 物件 | a_i | b_i | c_i | a_i^2 | b_i^2 | a_i b_i | a_i c_i | b_i c_i |
|---|---|---|---|---|---|---|---|---|
| A | −8 | −3 | −1.0 | 64 | 9 | 24 | 8.0 | 3.0 |
| B | −4 | 1 | −1.2 | 16 | 1 | −4 | 4.8 | −1.2 |
| C | −2 | −2 | 0.7 | 4 | 4 | 4 | −1.4 | −1.4 |
| D | 2 | 2 | 0.1 | 4 | 4 | 4 | 0.2 | 0.2 |
| E | 4 | −3 | 1.5 | 16 | 9 | −12 | 6.0 | −4.5 |
| F | 8 | 5 | −0.1 | 64 | 25 | 40 | −0.8 | −0.5 |
| 合計 | 0 | 0 | 0 | 168 | 52 | 56 | 16.8 | −4.4 |
合計の行から、S_{11} = \sum a_i^2 = 168、S_{22} = \sum b_i^2 = 52、S_{12} = \sum a_i b_i = 56、S_{1y} = \sum a_i c_i = 16.8、S_{2y} = \sum b_i c_i = -4.4 が得られる。偏差の合計がどれも0になっていることは、計算の検算に使える。
手順3（連立方程式を解く） 解くべき連立方程式は
である。公式の分母は S_{11} S_{22} - S_{12}^2 = 168 \times 52 - 56^2 = 8736 - 3136 = 5600 なので
手順4（切片を求める）
以上から、あてはめた式は
となる。たとえば広さ30㎡・駅徒歩5分の部屋なら、家賃は 5.0 + 0.2 \times 30 - 0.3 \times 5 = 9.5（万円）と予測される。6物件の予測値と残差は次のとおりである。
| 物件 | 家賃 y_i | 予測値 \hat{y}_i | 残差 e_i = y_i - \hat{y}_i |
|---|---|---|---|
| A | 7.0 | 7.3 | −0.3 |
| B | 6.8 | 6.9 | −0.1 |
| C | 8.7 | 8.2 | 0.5 |
| D | 8.1 | 7.8 | 0.3 |
| E | 9.5 | 9.7 | −0.2 |
| F | 7.9 | 8.1 | −0.2 |
残差の合計は -0.3 - 0.1 + 0.5 + 0.3 - 0.2 - 0.2 = 0 である。正規方程式の2つ目の式も
と確かめられる。広さだけの単回帰の場合は、単回帰の公式から \hat{\beta}_1 = \dfrac{S_{1y}}{S_{11}} = \dfrac{16.8}{168} = 0.1、\hat{\beta}_0 = 8.0 - 0.1 \times 30 = 5.0 となり、図1の直線 \hat{y} = 5.0 + 0.1x が得られる。
偏回帰係数の意味
広さだけの単回帰では傾きが0.1だったのに、駅徒歩を加えた重回帰では広さの係数が0.2と2倍になった。どちらが「広さの効果」なのだろうか。
重回帰の係数 \hat{\beta}_1 = 0.2 は、駅徒歩が同じ物件どうしで比べると、広さが1㎡大きいと家賃が平均0.2万円高いことを表す。同じように \hat{\beta}_2 = -0.3 は、広さが同じ物件どうしで比べると、駅徒歩が1分長いと家賃が平均0.3万円安いことを表す。ほかの説明変数を一定にしたときの、その変数だけの効果という意味で、偏回帰係数と呼ばれる。
図2は、駅徒歩を7分、11分、15分に固定したときの予測式 \hat{y} = 5.0 + 0.2 x_1 - 0.3 x_2 を、広さに対して描いたものである。徒歩時間ごとの直線はどれも傾き0.2で平行に並び、徒歩が4分長くなるごとに 0.3 \times 4 = 1.2 万円ずつ下にずれる。重回帰であてはめた平面を、駅徒歩の値ごとに切った断面がこれらの直線である。
単回帰の傾きが小さくなったのは、このデータでは広い物件ほど駅から遠い傾向があるからである（広さと駅徒歩の相関係数は0.60）。広さだけで比べると、「広いので高い」効果と「駅から遠いので安い」効果が混ざってしまう。式で確かめよう。広さだけの単回帰では、説明変数が x_1 だけなので連立方程式は S_{11} \hat{\beta}_1 = S_{1y} の1本になり、傾きは \dfrac{S_{1y}}{S_{11}} = \dfrac{16.8}{168} = 0.1 だった。一方、重回帰の連立方程式の1つ目の式 S_{11} \hat{\beta}_1 + S_{12} \hat{\beta}_2 = S_{1y} の両辺を S_{11} で割ると
となる。\dfrac{S_{12}}{S_{11}} = \dfrac{1}{3} は、駅徒歩を広さで単回帰したときの傾きで、「広さが1㎡大きい物件は、駅徒歩が平均 \dfrac{1}{3} 分長い」という関係を表す。つまり単回帰の傾き0.1は、広さそのものの効果0.2に、駅から遠くなることによる効果 (-0.3) \times \dfrac{1}{3} = -0.1 が上乗せされた値なのである。
偏回帰係数の値は、モデルにどの説明変数を入れるかによって変わる。また、偏回帰係数はデータの中での関係を表すものであり、それがそのまま因果関係（広さを変えれば家賃が変わる）を意味するとは限らない。モデルに入れていない要因（築年数など）が関係している可能性もある。
行列による表し方
説明変数が3つ以上になると、連立方程式を1つずつ書くのは大変になる。そこで行列を使ってまとめて書く。目的変数を縦に並べたベクトル Y、説明変数を並べた行列 X、係数のベクトル \beta を
とおき、誤差 \varepsilon_1, \dots, \varepsilon_n も同じように縦に並べたベクトルを \varepsilon とする。X は n 行 d+1 列の行列で、1列目がすべて1なのは、切片 \beta_0 に掛ける数だからである。X を計画行列（design matrix）という。すると n 個の式はまとめて
と書ける。残差平方和はベクトル Y - X\beta の長さの2乗 \| Y - X\beta \|^2 であり、これを最小にする \hat{\beta} が最小二乗推定量である。
正規方程式と最小二乗推定量
説明変数が2つの場合、残差平方和を最小にする条件は、次の3つの正規方程式だった（e_i = y_i - \hat{\beta}_0 - \hat{\beta}_1 x_{i1} - \hat{\beta}_2 x_{i2}）。
説明変数が d 個の場合も、残差平方和を \beta_0, \beta_1, \dots, \beta_d でそれぞれ偏微分して0とおくと、同じ形をした d + 1 個の式が得られる（ここでは e_i = y_i - \hat{\beta}_0 - \hat{\beta}_1 x_{i1} - \cdots - \hat{\beta}_d x_{id}）。
どの式も「残差 e_i に何かの数を掛けて合計すると0」という形をしている。掛ける数は、1つ目の式ではすべて1、2つ目以降は x_{i1}, \dots, x_{id} で、これはちょうど計画行列 X の各列である。そこで、残差を縦に並べたベクトルを e = Y - X\hat{\beta} とし、X の転置行列 X^\top（X の行と列を入れかえた行列）を e に掛けてみると
となり、各成分がちょうど正規方程式の左辺になっている。したがって、d + 1 個の正規方程式はまとめて X^\top e = 0 と書ける。e = Y - X\hat{\beta} を代入すると
となる。X^\top X は d+1 行 d+1 列の正方行列で、これが逆行列をもつとき、両辺に左から (X^\top X)^{-1} を掛けると次の公式が得られる。
家賃のデータでは、計画行列 X（各行が1つの物件で、1列目がすべて1、2列目が広さ、3列目が駅徒歩）と目的変数のベクトル Y は
で、これらから
である。X^\top X の成分はデータ数 n = 6、\sum x_{i1} = 180、\sum x_{i1}^2 = 5568、\sum x_{i1} x_{i2} = 1856 などの和、X^\top Y の成分は \sum y_i = 48、\sum x_{i1} y_i = 1456.8、\sum x_{i2} y_i = 475.6 である。逆行列を求めて掛けると \hat{\beta} = (5.0, \; 0.2, \; -0.3)^\top となり、計算例で連立方程式を解いて求めた \hat{\beta}_0 = 5.0、\hat{\beta}_1 = 0.2、\hat{\beta}_2 = -0.3 と一致する。
本当に最小になっているか
正規方程式は「偏微分が0」という条件から導いたので、本当に最小になっていることを確かめておこう。任意の \beta について Y - X\beta = e + X(\hat{\beta} - \beta) と分けると
となる（X^\top e = 0 なので真ん中の項が消える）。2つ目の項は0以上なので、\| Y - X\beta \|^2 は \beta = \hat{\beta} のとき最小値 \| e \|^2 をとる。さらに X^\top X が逆行列をもつとき（X の列が1次独立なとき）は、\beta \neq \hat{\beta} なら X(\hat{\beta} - \beta) \neq 0 となるので、最小にする \beta はただ1つに決まる。
説明変数の数がデータの数以上（d + 1 > n）のときや、ある説明変数がほかの説明変数の組み合わせで正確に表せるとき（たとえば「広さ（㎡）」と「広さ（坪）」を両方入れた場合）は、X^\top X が逆行列をもたない。このとき正規方程式を満たす \hat{\beta} は無数にあり、最小二乗推定量は1つに決まらない。正確に表せなくても、説明変数どうしの相関が非常に強い場合（多重共線性）は推定値が不安定になる。こうした場合の対処法として、リッジ回帰などの正則化がある。
射影としての見方
最小二乗法は、ベクトルの図形として眺めるとわかりやすい。データが n 個あれば、Y は n 次元空間の1つのベクトルである。一方、係数 \beta をいろいろに変えたときの
は、X の列ベクトルの組み合わせで作れるベクトル全体を動く。これは n 次元空間の中の平面のような部分空間で、\mathrm{Im}(X)（X の像）と書く。
最小二乗法は、\mathrm{Im}(X) の中で Y に最も近い点を探すことにあたる。図3のように、それは Y から \mathrm{Im}(X) に下ろした垂線の足、つまり Y の \mathrm{Im}(X) への正射影である。予測値 \hat{Y} = X\hat{\beta} がその点で、残差 e = Y - \hat{Y} は \mathrm{Im}(X) に垂直になる。正規方程式 X^\top e = 0（e が X のすべての列と直交する）は、まさにこの「垂直」を表している。
予測値は、\hat{\beta} の公式を代入して
と書ける。n 行 n 列の行列 P_X は Y を \mathrm{Im}(X) へ射影する行列で、Y に掛けると予測値 \hat{Y}（ハット）が得られることからハット行列とも呼ばれる。P_X は次の性質をもつ。
- P_X^2 = P_X（2回射影しても1回と同じ）。実際 P_X^2 = X(X^\top X)^{-1} X^\top X (X^\top X)^{-1} X^\top = X(X^\top X)^{-1} X^\top = P_X である
- P_X^\top = P_X（対称行列）
- P_X X = X（\mathrm{Im}(X) の中のベクトルは射影しても動かない）
残差は e = Y - P_X Y = (I - P_X) Y と書ける（I は単位行列）。I - P_X は \mathrm{Im}(X) に垂直な方向への射影行列で、同じく (I - P_X)^2 = I - P_X、(I - P_X)^\top = I - P_X を満たす。
また、X の1列目はすべて1なので、X^\top e = 0 の1行目から \sum e_i = 0 が出る。つまり、切片を含むモデルでは残差の合計が必ず0になる。家賃のデータの残差 -0.3, \; -0.1, \; 0.5, \; 0.3, \; -0.2, \; -0.2 の合計が0になったのも、このためである。
決定係数
変動の分解
あてはめたモデルがデータをどれくらいよく説明しているかを測りたい。そこで、目的変数の平均からのずれ y_i - \bar{y} を次のように2つに分ける。
図4は、家賃のデータの実測値 y を縦軸、予測値 \hat{y} を横軸にとった図で、物件 A についてこの分解を示している。点が対角線 y = \hat{y} に近いほど、予測がよく当たっている。説明変数がいくつあっても、この図は描ける。
両辺を2乗して全データで合計すると、交差項が消えて次の関係が成り立つ。
交差項は 2 \sum (\hat{y}_i - \bar{y}) e_i で、\sum \hat{y}_i e_i = \hat{Y}^\top e = \hat{\beta}^\top X^\top e = 0 と \sum \bar{y} e_i = \bar{y} \sum e_i = 0 から0になる。図形的には、すべて1のベクトルを \mathbf{1} として、ベクトル \hat{Y} - \bar{y}\mathbf{1} は \mathrm{Im}(X) の中にあり、e はそれに垂直なので、この関係は三平方の定理そのものである。
家賃のデータで確かめよう。平均は \bar{y} = 8.0、予測値は \hat{y} = 5.0 + 0.2 x_1 - 0.3 x_2 から計算したもので、各物件の3つのずれは次のとおりである。どの行でも、右の2列の和が y_i - \bar{y} になっている。
| 物件 | 家賃 y_i | 予測値 \hat{y}_i | y_i - \bar{y} | \hat{y}_i - \bar{y} | e_i = y_i - \hat{y}_i |
|---|---|---|---|---|---|
| A | 7.0 | 7.3 | −1.0 | −0.7 | −0.3 |
| B | 6.8 | 6.9 | −1.2 | −1.1 | −0.1 |
| C | 8.7 | 8.2 | 0.7 | 0.2 | 0.5 |
| D | 8.1 | 7.8 | 0.1 | −0.2 | 0.3 |
| E | 9.5 | 9.7 | 1.5 | 1.7 | −0.2 |
| F | 7.9 | 8.1 | −0.1 | 0.1 | −0.2 |
それぞれの列の2乗の合計を求めると
となり、5.20 = 4.68 + 0.52 が成り立っている。なお、回帰変動は偏差の積和からも計算できる。交差項が0であること（\sum (\hat{y}_i - \bar{y}) e_i = 0）を使うと
となり、ここに \hat{y}_i - \bar{y} = \hat{\beta}_1 (x_{i1} - \bar{x}_1) + \hat{\beta}_2 (x_{i2} - \bar{x}_2) を代入すると、計算例で求めた S_{1y} = 16.8、S_{2y} = -4.4 を使って
が得られる。
変動の分解は、\sum e_i = 0 を使って導いた。切片を含まないモデル（X にすべて1の列がないモデル）では残差の合計が0になるとは限らず、この分解は一般には成り立たない。
決定係数 R²
総変動のうち回帰変動が占める割合を決定係数（coefficient of determination）といい、R^2 で表す。
変動の分解から 0 \leq R^2 \leq 1 であり、1に近いほどデータへのあてはまりがよい。家賃のデータでは R^2 = \dfrac{4.68}{5.20} = 0.90 で、家賃のばらつきの90%を広さと駅徒歩で説明できている。広さだけの単回帰では、残差 -0.2, -0.8, 0.9, -0.1, 1.1, -0.9 から残差変動が3.52で、R^2 = 1 - \dfrac{3.52}{5.20} = 0.32 だった。駅徒歩を加えたことで、説明力が大きく上がったことがわかる。
R^2 は、実測値 y_i と予測値 \hat{y}_i の相関係数の2乗に等しいことも知られている。この相関係数 R を重相関係数という。単回帰の場合、R^2 は x と y の相関係数の2乗になる。
自由度調整済み決定係数
R^2 には弱点がある。説明変数を増やすと、R^2 は決して下がらない。新しい変数の係数を0にすれば元のモデルと同じになるので、変数を増やしたモデルの残差変動の最小値は、元のモデルの残差変動以下になるからである。家賃とまったく関係のない変数（たとえば部屋番号）を加えても、R^2 はふつう少し上がってしまう。
そこで、説明変数の数を考慮した自由度調整済み決定係数（adjusted R^2）がよく使われる。
分子は後で示すように誤差分散 \sigma^2 の不偏推定量、分母は y の分散の不偏推定量である。説明変数を増やすと残差変動は減るが、割る数 n - d - 1 も小さくなる。そのため、あまり役に立たない変数を加えると R^{*2} はかえって下がる。家賃のデータでは
である。実際のデータで R^2 と R^{*2} が逆向きに動く例は、後の具体例で紹介する。
最小二乗推定量の性質
ここからは、データが本当に重回帰モデル Y = X\beta + \varepsilon から生まれていると仮定して、最小二乗推定量の性質を調べる。\hat{\beta} はデータ Y から計算されるので、誤差 \varepsilon の出方によって値が変わる確率変数である。その平均と分散を求めよう。以下では説明変数の値（X）は固定されたものとし、期待値は誤差についてとる。
不偏性
\hat{\beta} の式に Y = X\beta + \varepsilon を代入すると
となる。\hat{\beta} は、真の値 \beta に誤差の重みつき和を加えたものである。E[\varepsilon] = 0 なので
となり、最小二乗推定量は不偏推定量である。
分散共分散行列
\hat{\beta} - \beta = (X^\top X)^{-1} X^\top \varepsilon なので、\hat{\beta} の分散共分散行列は
となる。E[\varepsilon \varepsilon^\top] の (i, j) 成分は E[\varepsilon_i \varepsilon_j] で、i = j なら \sigma^2、i \neq j なら誤差が無相関なので0である。よって E[\varepsilon \varepsilon^\top] = \sigma^2 I を使った。
対角成分が各係数の推定量の分散、対角以外の成分が推定量どうしの共分散である。家賃のデータでは
なので、V[\hat{\beta}_1] = 0.009286 \, \sigma^2、V[\hat{\beta}_2] = 0.03 \, \sigma^2、\mathrm{Cov}[\hat{\beta}_1, \hat{\beta}_2] = -0.01 \, \sigma^2 である。
この行列の右下の2行2列の部分（偏回帰係数に対応する部分）は、偏回帰係数を求める連立方程式の係数行列の逆行列
に一致している（これは一般に成り立つ）。したがって、x_1 と x_2 の相関係数を r_{12} = \dfrac{S_{12}}{\sqrt{S_{11} S_{22}}} とすると
と書ける。家賃のデータでは r_{12}^2 = \dfrac{56^2}{168 \times 52} = 0.359 なので、\dfrac{1}{168 \times (1 - 0.359)} = 0.009286 となり、(X^\top X)^{-1} の対応する成分 0.009286 と一致する。説明変数どうしの相関が強いほど（r_{12}^2 が1に近いほど）、偏回帰係数の推定量のばらつきは大きくなる。これが、多重共線性があると推定が不安定になる理由である。
図5は、6物件の広さと駅徒歩はそのままにして、真の値を \beta = (5, \; 0.2, \; -0.3)^\top、\sigma = 0.4 として家賃のデータを10000回つくり直し、そのたびに \hat{\beta}_1 を計算した結果である。\hat{\beta}_1 は真の値0.2のまわりに分布し、10000回の平均は0.2005、標準偏差は0.0383であった。理論上の標準偏差 0.4 \times \sqrt{0.009286} = 0.0385 とよく一致している。
誤差が正規分布に従うとき、\hat{\beta} は誤差の1次式なので、\hat{\beta} 自身も多変量正規分布 N\bigl(\beta, \; \sigma^2 (X^\top X)^{-1}\bigr) に従う。図5の曲線はその1成分の分布である。実際の分析では \sigma^2 は未知なので、後で述べる不偏推定量 s^2 で置き換えて推定量のばらつき（標準誤差）を見積もる。これが回帰係数の検定や信頼区間の基礎になる。
ちなみに、同じ10000回のデータで広さだけの単回帰の傾きを計算すると、その平均は0.1002となり、真の値0.2から大きくずれる。このデータでは広い物件ほど駅から遠い傾向がある（広さが1㎡大きいと駅徒歩が平均 \dfrac{1}{3} 分長い）ため、広さだけの単回帰の傾きには駅徒歩の効果 (-0.3) \times \dfrac{1}{3} = -0.1 が混ざり、平均すると 0.2 - 0.1 = 0.1 になってしまうのである。一般に、目的変数に影響し、しかもほかの説明変数と相関のある変数をモデルから落とすと、残した変数の係数の推定量に偏りが生じる（欠落変数バイアス）。
最小分散性
では、\hat{\beta} の分散 \sigma^2 (X^\top X)^{-1} はどれくらい小さいのだろうか。誤差が正規分布に従うとき、対数尤度は
である。\beta を含むのは第2項だけなので、\beta の最尤推定量は残差平方和を最小にする \beta、つまり最小二乗推定量と一致する。さらに \ell を \beta で微分すると
となるので、\beta についてのフィッシャー情報行列は \dfrac{X^\top X}{\sigma^2} であり、クラメール・ラオの不等式が与える分散の下限はその逆行列 \sigma^2 (X^\top X)^{-1} である。これは \hat{\beta} の分散共分散行列そのものなので、最小二乗推定量は下限を達成している。つまり、どんな不偏推定量 \tilde{\beta} をもってきても、各係数の分散は \hat{\beta} の分散より小さくならない。最小二乗推定量は最小分散不偏推定量である。
誤差が正規分布に従うと仮定しない場合でも、誤差の平均が0、分散が等しく、互いに無相関であれば、最小二乗推定量は「Y の1次式で表される不偏推定量（線形不偏推定量）」の中で分散が最小になる。これをガウス・マルコフの定理といい、最小二乗推定量は最良線形不偏推定量（BLUE：best linear unbiased estimator）であるという。ガウスは1821年に発表した最小二乗法の理論の中でこの定理の一つの形を示しており、のちにマルコフ連鎖に名を残すアンドレイ・マルコフ（Andrey Markov）が仮定を弱めた形で述べたことから、2人の名前がついている。
誤差分散の不偏推定量
誤差の分散 \sigma^2 は残差から推定する。誤差 \varepsilon_i そのものは観測できないので、代わりに残差 e_i を使う。残差平方和を n で割ればよさそうだが、実際には n - d - 1 で割るのが正しい。
これを示そう。P_X X = X (X^\top X)^{-1} X^\top X = X より (I - P_X) X = X - X = O なので、残差は
と誤差だけで表せる。I - P_X は対称で、P_X^2 = P_X より (I - P_X)^2 = I - 2P_X + P_X^2 = I - P_X を満たすので、残差平方和は e^\top e = \varepsilon^\top (I - P_X) \varepsilon となる。一般に行列 A について E[\varepsilon^\top A \varepsilon] = \sum_{i} \sum_{j} A_{ij} E[\varepsilon_i \varepsilon_j] = \sigma^2 \sum_{i} A_{ii} = \sigma^2 \, \mathrm{tr}(A)（\mathrm{tr} は対角成分の和）なので
である。ここでトレースの性質 \mathrm{tr}(AB) = \mathrm{tr}(BA) を使うと
となり（I_{d+1} は d+1 次の単位行列）、E\bigl[\sum e_i^2\bigr] = (n - d - 1) \sigma^2 が得られる。したがって s^2 は \sigma^2 の不偏推定量である。
n - d - 1 で割る理由は、標本分散を n - 1 で割るのと同じである。残差 e_1, \dots, e_n は正規方程式 X^\top e = 0 という d + 1 個の条件を満たすので、自由に動けるのは n - d - 1 個分しかない。推定した係数の個数だけ自由度が減るのである。家賃のデータでは、残差変動0.52から s^2 = \dfrac{0.52}{6 - 2 - 1} = 0.173、s = 0.416（万円）である。なお、誤差に正規分布を仮定したときの \sigma^2 の最尤推定量は残差平方和を n で割ったもので、これは \sigma^2 を平均的に小さめに見積もる。
多項式回帰：曲線をあてはめる
重回帰分析は直線や平面しかあてはめられないように見えるが、説明変数の選び方を工夫すれば曲線もあてはめられる。たとえば説明変数が x の1つだけでも、x^2 を2つ目の説明変数とみなして
とすれば、x_1 = x、x_2 = x^2 とした重回帰モデルそのものである。一般に、関数 \varphi_1(x), \dots, \varphi_d(x) を用意して
とするモデルも、x_k = \varphi_k(x) と置き換えれば重回帰モデルになる。\varphi_k を基底関数といい、\varphi_k(x) = x^k とすれば多項式回帰、\varphi_k(x) = \cos(kx) などとすれば三角関数による回帰になる。描く曲線は直線でなくても、係数 \beta について1次式であれば、公式 \hat{\beta} = (X^\top X)^{-1} X^\top Y がそのまま使える。「線形回帰」の「線形」は、x についてではなく、係数について1次式という意味である。目的変数や説明変数に対数変換をほどこしてから重回帰を行うのも、同じ考え方である。
例として、ある自動車の速度と燃費を測った次のデータを考える。
| 速度 x（km/h） | 20 | 30 | 40 | 50 | 60 | 70 | 80 | 90 | 100 |
|---|---|---|---|---|---|---|---|---|---|
| 燃費 y（km/L） | 13.6 | 16.2 | 17.7 | 18.5 | 19.0 | 18.5 | 17.8 | 16.5 | 13.4 |
直線 y = \beta_0 + \beta_1 x をあてはめると \hat{y} = 16.770 + 0.0005x とほぼ水平になり、R^2 は0.000である。速度が低すぎても高すぎても燃費が悪くなる山なりの関係は、直線ではまったくとらえられない。2次式をあてはめると
となり、データによく合う（図6）。燃費が最大になる速度は、2次関数の頂点 x = -\dfrac{\hat{\beta}_1}{2\hat{\beta}_2} から約60 km/h と求められる。
多項式の次数を上げると、R^2 はいくらでも1に近づく。9個のデータなら8次式で全点をちょうど通る曲線が作れ、R^2 = 1 となる。しかしそのような曲線はデータの細かなばらつきまで追いかけたもので、新しいデータの予測には役に立たない（過学習）。変数の選び方には、自由度調整済み決定係数や AIC などの基準、あるいは正則化が使われる。
具体例：化学プラントのデータ
実際のデータで重回帰分析を行ってみよう。使うのは、アンモニアを酸化して硝酸をつくる化学プラントの21日分の運転データで、ブラウンリー（K. A. Brownlee）が1960年の著書 Statistical Theory and Methodology in Science and Engineering で紹介したものである。多くの手法の試験台として使われてきた有名なデータで、「重回帰分析のモルモット」（The Guinea Pig of Multiple Regression）という題の論文（Dodge, 1996）まである。統計ソフト R や Python の statsmodels にも stackloss という名前で収録されている。変数は次の4つである。
- 空気流量：プラントの稼働の度合い
- 冷却水温：吸収塔のコイルを流れる冷却水の温度
- 酸濃度：循環する酸の濃度（%）から50を引いて10倍した値（89なら58.9%）
- スタックロス（目的変数）：プラントに入ったアンモニアのうち、吸収されずに逃げた割合（%）の10倍。値が大きいほどプラントの効率が悪い
| 日 | 空気 流量 | 冷却 水温 | 酸 濃度 | スタック ロス | 日 | 空気 流量 | 冷却 水温 | 酸 濃度 | スタック ロス |
|---|---|---|---|---|---|---|---|---|---|
| 1 | 80 | 27 | 89 | 42 | 12 | 58 | 17 | 88 | 13 |
| 2 | 80 | 27 | 88 | 37 | 13 | 58 | 18 | 82 | 11 |
| 3 | 75 | 25 | 90 | 37 | 14 | 58 | 19 | 93 | 12 |
| 4 | 62 | 24 | 87 | 28 | 15 | 50 | 18 | 89 | 8 |
| 5 | 62 | 22 | 87 | 18 | 16 | 50 | 18 | 86 | 7 |
| 6 | 62 | 23 | 87 | 18 | 17 | 50 | 19 | 72 | 8 |
| 7 | 62 | 24 | 93 | 19 | 18 | 50 | 19 | 79 | 8 |
| 8 | 62 | 24 | 93 | 20 | 19 | 50 | 20 | 80 | 9 |
| 9 | 58 | 23 | 87 | 15 | 20 | 56 | 20 | 82 | 15 |
| 10 | 58 | 18 | 80 | 14 | 21 | 70 | 20 | 91 | 15 |
| 11 | 58 | 18 | 89 | 14 |
空気流量だけを使う単回帰から始めて、冷却水温、酸濃度と説明変数を1つずつ加えた3つのモデルをあてはめると、次のようになった（s は誤差の標準偏差の推定値 \sqrt{s^2}）。
| モデル | 切片 | 空気流量 | 冷却水温 | 酸濃度 | R^2 | R^{*2} | s |
|---|---|---|---|---|---|---|---|
| (1) 空気流量 | −44.132 | 1.020 | − | − | 0.8458 | 0.8377 | 4.098 |
| (2) (1)＋冷却水温 | −50.359 | 0.671 | 1.295 | − | 0.9088 | 0.8986 | 3.239 |
| (3) (2)＋酸濃度 | −39.920 | 0.716 | 1.295 | −0.152 | 0.9136 | 0.8983 | 3.243 |
この結果から、次の3つのことが読み取れる。
(a) 冷却水温は説明に役立つ 冷却水温を加えると R^2 は0.846から0.909に上がり、誤差の標準偏差の推定値 s も4.10から3.24に下がった。
(b) 空気流量の係数は単回帰より小さくなる 空気流量の係数は、単回帰の1.020から (2) では0.671に下がった。空気流量と冷却水温の相関係数は0.78と高く、空気流量の多い日は冷却水温も高い傾向がある。単回帰の傾き1.020には冷却水温の効果も混ざっていたのであり、家賃の例で、広さだけの単回帰の傾き0.1が、駅徒歩を加えると0.2になったのと同じ現象である。(3) のモデルでは「冷却水温と酸濃度が同じ日どうしで比べると、空気流量が1多いとスタックロスが平均0.716多い」と解釈する。
(c) 酸濃度はほとんど役立たない 酸濃度を加えると R^2 は0.9088から0.9136へわずかに上がったが、自由度調整済み決定係数 R^{*2} は0.8986から0.8983へと逆に下がった。s もわずかに大きくなっている。説明変数を増やせば R^2 は必ず上がるが、それに見合うほど残差が減っていないということである。
酸濃度の係数 −0.152 が、偶然0からずれただけなのか、本当に効果があるのかを判断するには、回帰係数の検定（t分布やF分布を使う）が必要になる。これは別の記事で扱う。
練習問題
[2] 4人世帯で床面積100㎡の家の電気代を予測せよ。
[3] 世帯人数だけで単回帰したときの傾きを求め、[1] の \hat{\beta}_1 と異なる理由を説明せよ。
[1] 連立方程式は
である。下の式から上の式を引くと 30 \hat{\beta}_2 = 150 より \hat{\beta}_2 = 5、上の式に代入して \hat{\beta}_1 = 20 - 5 = 15 である（公式を使っても、分母 10 \times 40 - 10^2 = 300 から同じ値になる）。切片は
[2] 床面積100㎡は x_2 = 10 なので
つまり約12500円と予測される。
[3] 単回帰の傾きは \dfrac{S_{1y}}{S_{11}} = \dfrac{200}{10} = 20 で、偏回帰係数15より大きい。S_{12} = 10 > 0 なので、人数の多い世帯ほど床面積が広い傾向がある。\dfrac{S_{12}}{S_{11}} = 1 より、人数が1人多い世帯は床面積が平均10㎡広く、その分の電気代 5 \times 1 = 5 も単回帰の傾きに含まれてしまう。実際 \hat{\beta}_1 + \hat{\beta}_2 \times \dfrac{S_{12}}{S_{11}} = 15 + 5 \times 1 = 20 である。
[1] 各モデルの決定係数を求めよ。
[2] 各モデルの自由度調整済み決定係数を求め、どちらのモデルを選ぶのがよいか答えよ。
[3] モデル A について、誤差分散 \sigma^2 の不偏推定値を求めよ。
[1]
[2] モデル A は n - d - 1 = 20 - 2 - 1 = 17、モデル B は 20 - 5 - 1 = 14 なので
決定係数はモデル B のほうが大きいが、自由度調整済み決定係数はモデル A のほうが大きい。B で増やした3つの変数は、残差変動を100から90に減らしただけで、変数を増やした分に見合わない。モデル A を選ぶのがよい。
[3]
[1] 計画行列 X について X^\top X と X^\top Y を求め、\hat{\beta} = (X^\top X)^{-1} X^\top Y を計算せよ。
[2] 残差 e_1, \dots, e_4 を求め、\sum e_i = 0 と \sum x_i e_i = 0 が成り立つことを確かめよ。
[3] 誤差分散の不偏推定値 s^2 を求めよ。また、\mathrm{Cov}[\hat{\beta}] = \sigma^2 (X^\top X)^{-1} の \sigma^2 を s^2 で置き換えて、\hat{\beta}_1 の分散を推定せよ。
[1] 計画行列は1列目がすべて1、2列目が x なので
である。4 \times 30 - 10 \times 10 = 20 なので
あてはめた直線は \hat{y} = 0.5 + 1.4x である。
[2] 予測値は 1.9, \; 3.3, \; 4.7, \; 6.1 なので、残差は 0.1, \; -0.3, \; 0.3, \; -0.1 である。
[3] 残差変動は 0.01 + 0.09 + 0.09 + 0.01 = 0.2、n - d - 1 = 4 - 1 - 1 = 2 なので s^2 = \dfrac{0.2}{2} = 0.1 である。(X^\top X)^{-1} の (2, 2) 成分は \dfrac{4}{20} = 0.2 なので
となる。その平方根 \sqrt{0.02} = 0.141 が \hat{\beta}_1 の標準誤差である。
まとめ
| 項目 | 内容 |
|---|---|
| モデル | y_i = \beta_0 + \beta_1 x_{i1} + \cdots + \beta_d x_{id} + \varepsilon_i（\varepsilon_i は互いに独立に N(0, \sigma^2)） |
| 行列による表記 | Y = X\beta + \varepsilon（X は1列目がすべて1の n 行 d+1 列の計画行列） |
| 最小二乗推定量 | 残差平方和 \| Y - X\beta \|^2 を最小にする \hat{\beta} = (X^\top X)^{-1} X^\top Y |
| 正規方程式 | X^\top X \hat{\beta} = X^\top Y（残差と X の各列が直交：X^\top e = 0） |
| 説明変数が2つの場合 | S_{11} \hat{\beta}_1 + S_{12} \hat{\beta}_2 = S_{1y}、S_{12} \hat{\beta}_1 + S_{22} \hat{\beta}_2 = S_{2y}、\hat{\beta}_0 = \bar{y} - \hat{\beta}_1 \bar{x}_1 - \hat{\beta}_2 \bar{x}_2 |
| 単回帰の場合 | \hat{\beta}_1 = \dfrac{S_{1y}}{S_{11}}（共分散 ÷ 分散）、\hat{\beta}_0 = \bar{y} - \hat{\beta}_1 \bar{x}_1 |
| 偏回帰係数 | ほかの説明変数を一定にしたときの効果。説明変数どうしに相関があると、単回帰の傾きとは異なる |
| 射影 | \hat{Y} = P_X Y、P_X = X(X^\top X)^{-1} X^\top（ハット行列）。残差 e = (I - P_X) Y は \mathrm{Im}(X) に垂直 |
| 変動の分解 | 総変動 ＝ 回帰変動 ＋ 残差変動（切片を含むモデル） |
| 決定係数 | R^2 = 1 - \dfrac{\text{残差変動}}{\text{総変動}}（説明変数を増やすと下がらない） |
| 自由度調整済み決定係数 | R^{*2} = 1 - (1 - R^2) \times \dfrac{n - 1}{n - d - 1} |
| 推定量の性質 | E[\hat{\beta}] = \beta、\mathrm{Cov}[\hat{\beta}] = \sigma^2 (X^\top X)^{-1}。誤差が正規分布なら最尤推定量と一致し、最小分散不偏推定量。正規分布を仮定しなくても最良線形不偏推定量（ガウス・マルコフの定理） |
| \sigma^2 の推定 | s^2 = \dfrac{1}{n - d - 1} \sum e_i^2 が不偏推定量 |
| 曲線のあてはめ | x^2 や \cos(kx) などを説明変数とすれば、係数について1次式のまま曲線をあてはめられる（多項式回帰） |
Python実装
切片つきの重回帰を最小二乗法で推定する関数を作り、本文の計算を確かめる。正規方程式は np.linalg.solve で解いている（逆行列を求めてから掛けるより、数値計算として安定している）。
import numpy as np
def fit_ols(X, y):
"""説明変数の行列 X（n×d）と目的変数 y で、切片つき重回帰を推定する"""
n, d = X.shape
Z = np.column_stack([np.ones(n), X]) # 計画行列（1列目はすべて1）
beta = np.linalg.solve(Z.T @ Z, Z.T @ y) # 正規方程式を解く
y_hat = Z @ beta # 予測値
e = y - y_hat # 残差
rss = e @ e # 残差変動
tss = np.sum((y - y.mean())**2) # 総変動
s2 = rss / (n - d - 1) # σ² の不偏推定量
return {
"beta": beta, "y_hat": y_hat, "resid": e, "s2": s2,
# 決定係数、自由度調整済み決定係数、推定量の標準誤差
"r2": 1 - rss / tss,
"adj_r2": 1 - (rss / (n - d - 1)) / (tss / (n - 1)),
"se": np.sqrt(s2 * np.diag(np.linalg.inv(Z.T @ Z))),
}from multiple_regression import *
# ========== 家賃のデータ ==========
area = np.array([22, 26, 28, 32, 34, 38]) # 広さ（㎡）
walk = np.array([7, 11, 8, 12, 7, 15]) # 駅徒歩（分）
rent = np.array([7.0, 6.8, 8.7, 8.1, 9.5, 7.9]) # 家賃（万円）
res = fit_ols(np.column_stack([area, walk]), rent)
print("【家賃のデータ】")
print("回帰係数（β0, β1, β2）:", res["beta"].round(4))
print("予測値:", res["y_hat"].round(4))
print("残差:", res["resid"].round(4))
print("残差の和:", round(res["resid"].sum(), 10))
print(f"R² = {res['r2']:.4f}, 自由度調整済み R² = {res['adj_r2']:.4f}")
print(f"s² = {res['s2']:.4f}")
simple = fit_ols(area[:, None], rent)
print("広さだけの単回帰:", simple["beta"].round(4))
print(f" R² = {simple['r2']:.4f}")
# ハット行列 P = Z(ZᵀZ)⁻¹Zᵀ の性質
Z = np.column_stack([np.ones(6), area, walk])
P = Z @ np.linalg.inv(Z.T @ Z) @ Z.T
print("P² = P:", np.allclose(P @ P, P), " Pᵀ = P:", np.allclose(P.T, P))
print("tr(P) =", round(np.trace(P), 6))
# ========== 不偏性と分散の確認 ==========
# 同じ6物件で家賃のデータを10000回つくり直す
rng = np.random.default_rng(1)
beta_true, sigma = np.array([5.0, 0.2, -0.3]), 0.4
Y_sim = Z @ beta_true + rng.normal(0, sigma, (10000, 6)) # 各行が1回分
B = np.linalg.solve(Z.T @ Z, Z.T @ Y_sim.T).T # 各行が1回分の推定値
V = sigma**2 * np.linalg.inv(Z.T @ Z) # 理論上の分散共分散行列
print("\n【シミュレーション（σ = 0.4, 10000回）】")
print("推定値の平均:", B.mean(axis=0).round(4), "（真の値 [5, 0.2, -0.3]）")
print(f"β1 の推定値の分散: {B[:, 1].var():.6f}（理論値 {V[1, 1]:.6f}）")
print(f"β2 の推定値の分散: {B[:, 2].var():.6f}（理論値 {V[2, 2]:.6f}）")
c12 = np.cov(B[:, 1], B[:, 2])[0, 1]
print(f"β1 と β2 の推定値の共分散: {c12:.6f}（理論値 {V[1, 2]:.6f}）")
# ========== 多項式回帰：速度と燃費 ==========
speed = np.arange(20, 101, 10)
fuel = np.array([13.6, 16.2, 17.7, 18.5, 19.0, 18.5, 17.8, 16.5, 13.4])
lin = fit_ols(speed[:, None], fuel)
quad = fit_ols(np.column_stack([speed, speed**2]), fuel) # x と x² を使う
print("\n【速度と燃費】")
b0, b1 = lin["beta"]
print(f"直線: ŷ = {b0:.3f} + {b1:.4f}x, R² = {lin['r2']:.4f}")
b0, b1, b2 = quad["beta"]
print(f"2次式: ŷ = {b0:.3f} + {b1:.4f}x - {-b2:.5f}x²")
print(f" R² = {quad['r2']:.4f}")
print(f"燃費が最大になる速度: {-b1 / (2 * b2):.1f} km/h")
# ========== スタックロスのデータ（Brownlee, 1960） ==========
air = np.array([80, 80, 75, 62, 62, 62, 62, 62, 58, 58, 58,
58, 58, 58, 50, 50, 50, 50, 50, 56, 70])
water = np.array([27, 27, 25, 24, 22, 23, 24, 24, 23, 18, 18,
17, 18, 19, 18, 18, 19, 19, 20, 20, 20])
acid = np.array([89, 88, 90, 87, 87, 87, 93, 93, 87, 80, 89,
88, 82, 93, 89, 86, 72, 79, 80, 82, 91])
loss = np.array([42, 37, 37, 28, 18, 18, 19, 20, 15, 14, 14,
13, 11, 12, 8, 7, 8, 8, 9, 15, 15])
print("\n【スタックロスのデータ】")
models = [("(1) 空気流量", [air]),
("(2) ＋冷却水温", [air, water]),
("(3) ＋酸濃度", [air, water, acid])]
for name, cols in models:
r = fit_ols(np.column_stack(cols), loss)
print(f"{name}: 係数 {r['beta'].round(4)}")
print(f" R² = {r['r2']:.4f}, 自由度調整済み R² = {r['adj_r2']:.4f},"
f" s = {np.sqrt(r['s2']):.3f}")
# statsmodels でも同じ結果になることを確認
import statsmodels.api as sm
Xs = sm.add_constant(np.column_stack([air, water, acid]))
model = sm.OLS(loss, Xs).fit()
print("statsmodels:", model.params.round(4))
r2, adj = model.rsquared, model.rsquared_adj
print(f" R² = {r2:.4f}, 自由度調整済み R² = {adj:.4f}")家賃のデータでは、係数 (5.0, \; 0.2, \; -0.3)、R^2 = 0.9、自由度調整済み R^2 = 0.8333、s^2 = 0.1733 と本文の値が再現され、ハット行列について P^2 = P、P^\top = P、\mathrm{tr}(P) = 3 = d + 1 も確かめられる。シミュレーションでは、\hat{\beta} の平均が真の値にほぼ一致し（不偏性）、分散と共分散も \sigma^2 (X^\top X)^{-1} の値に近い。燃費のデータでは2次式の R^2 が0.9875で、燃費が最大になる速度は約60.1 km/h である。スタックロスのデータでは、自作の関数と statsmodels の結果が一致している。