Last active
March 6, 2019 22:59
-
-
Save shotahorii/9b1e7533e386453bb7894fc57bf6b9eb to your computer and use it in GitHub Desktop.
Ridge_Multicollinearity
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| { | |
| "cells": [ | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "### Why Ridge Regularization is effective for multicollinearity\n", | |
| "1. [リッジ回帰による多重共線性の問題回避について](http://jojoshin.hatenablog.com/entry/2017/05/01/011132)\n", | |
| "2. [リッジ回帰(ridge)で解決できる多重共線性問題と注意すべき点](https://www.bananarian.net/entry/multico-ridge)\n", | |
| "3. [Rでスパースモデリング：Elastic Net回帰についてまとめてみる](http://tekenuko.hatenablog.com/entry/2017/11/18/214317)\n", | |
| " \n", | |
| "[ここ](https://gist.github.com/shotahorii/64e772fa6d5d767a1fc1bcc87ee40353)で、多重共線性の問題について書いた。具体的には以下の2つのケースに分けられる。 \n", | |
| "1. 変数間に完全な線形関係が存在するケース。この場合、最小二乗法で解を求めることができず、推定量を一意に算出できない。\n", | |
| "2. 変数間に完全な線形関係はないが強い相関があるケース。この場合、推定量の算出は可能であるが、推定量の分散が大きくなり、推定量の値が安定しない。（この推定量の正しさが信頼できない）\n", | |
| "\n", | |
| "[上の記事](https://gist.github.com/shotahorii/64e772fa6d5d767a1fc1bcc87ee40353)では、1.のケースについて二変数のモデルを例に正規方程式を展開して説明したが、ここでは一般化して多変数のモデルにおいて変数間に完全な線形関係がある場合に解が求まらない理由を行列式で見ていく。 \n", | |
| " \n", | |
| "以下にでてくる最小二乗法の計算については[ここ](https://gist.github.com/shotahorii/037dc487ab3c70cc45af27fc9d03c9fb#file-leastsquares-ipynb)参照。\n", | |
| "\n", | |
| "#### 1. 変数間に完全な線形関係が存在するケース" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 51, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/latex": [ | |
| "最小二乗法の解は以下のように表される。\n", | |
| "<br><br>\n", | |
| "$$\\beta = (X^TX)^{-1}X^Ty$$\n", | |
| "<br><br>\n", | |
| "ここで、変数間に完全な線形関係がある場合$$(X^TX)^{-1}$$の計算が不可能となる。<br>\n", | |
| "以下に3変数、3サンプルのデータを例に見ていく。<br>\n", | |
| "計画行列Xにおいて、変数$$x_1,x_2$$が次のように完全に線形関係であると想定する。($$x_2 = 3x_1$$)ここで、\n", | |
| "$$x_3$$は特に関連のない変数とする。\n", | |
| "<br><br>\n", | |
| "$$X = \\left(\n", | |
| " \\begin{array}{cccc}\n", | |
| " 1 & 3 & 3 \\\\\n", | |
| " 2 & 6 & 12 \\\\\n", | |
| " 5 & 15 & 5\n", | |
| " \\end{array}\n", | |
| " \\right)$$\n", | |
| "<br><br>\n", | |
| "$$X^T = \\left(\n", | |
| " \\begin{array}{cccc}\n", | |
| " 1 & 2 & 5 \\\\\n", | |
| " 3 & 6 & 15 \\\\\n", | |
| " 3 & 12 & 5\n", | |
| " \\end{array}\n", | |
| " \\right)$$\n", | |
| "<br><br>\n", | |
| "$$X^TX = \\left(\n", | |
| " \\begin{array}{cccc}\n", | |
| " 1 & 2 & 5 \\\\\n", | |
| " 3 & 6 & 15 \\\\\n", | |
| " 3 & 12 & 5\n", | |
| " \\end{array}\n", | |
| " \\right)\n", | |
| " \\left(\n", | |
| " \\begin{array}{cccc}\n", | |
| " 1 & 3 & 3 \\\\\n", | |
| " 2 & 6 & 12 \\\\\n", | |
| " 5 & 15 & 5\n", | |
| " \\end{array}\n", | |
| " \\right)\n", | |
| "= \\left(\n", | |
| " \\begin{array}{cccc}\n", | |
| " 30 & 90 & 52 \\\\\n", | |
| " 90 & 270 & 156 \\\\\n", | |
| " 52 & 156 & 178\n", | |
| " \\end{array}\n", | |
| " \\right) \n", | |
| "$$\n", | |
| "<br><br>\n", | |
| "ここで、3×3の行列Aの逆行列は以下のように求められる。\n", | |
| "<br><br>\n", | |
| "$$A=\\left(\n", | |
| " \\begin{array}{cccc}\n", | |
| " a_{11} & a_{12} & a_{13} \\\\\n", | |
| " a_{21} & a_{22} & a_{23} \\\\\n", | |
| " a_{31} & a_{32} & a_{33}\n", | |
| " \\end{array}\n", | |
| " \\right)$$\n", | |
| "<br>\n", | |
| "とおくと、\n", | |
| "<br>\n", | |
| "$$A^{-1}=\\frac{1}{detA} \\left(\n", | |
| " \\begin{array}{cccc}\n", | |
| " a_{22}a_{33}-a_{23}a_{32} & a_{13}a_{32}-a_{12}a_{33} & a_{12}a_{23}-a_{13}a_{22} \\\\\n", | |
| " a_{23}a_{31}-a_{21}a_{33} & a_{11}a_{33}-a_{13}a_{31} & a_{13}a_{21}-a_{11}a_{23} \\\\\n", | |
| " a_{21}a_{32}-a_{22}a_{31} & a_{12}a_{31}-a_{11}a_{32} & a_{11}a_{22}-a_{12}a_{21}\n", | |
| " \\end{array}\n", | |
| " \\right)$$\n", | |
| "<br>\n", | |
| "ここで、\n", | |
| "<br>\n", | |
| "$$detA = a_{11}a_{22}a_{33} + a_{21}a_{32}a_{13} + a_{31}a_{12}a_{23}\n", | |
| " - a_{11}a_{32}a_{23} - a_{31}a_{22}a_{13} - a_{21}a_{12}a_{33}$$\n", | |
| "<br><br>\n", | |
| "これを元に、先の計画行列のdetAを求めてみる。\n", | |
| "<br><br>\n", | |
| "$$detA = 30\\cdot 270\\cdot 178 + 90\\cdot 156\\cdot 52 + 52\\cdot 90\\cdot 156\n", | |
| " - 30\\cdot 156\\cdot 156 - 52\\cdot 270\\cdot 52 - 90\\cdot 90\\cdot 178$$\n", | |
| "$$= 1441800 + 730080 + 730080 - 730080 - 730080 - 1441800 = 0$$" | |
| ], | |
| "text/plain": [ | |
| "<IPython.core.display.Latex object>" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "%%latex\n", | |
| "最小二乗法の解は以下のように表される。\n", | |
| "<br><br>\n", | |
| "$$\\beta = (X^TX)^{-1}X^Ty$$\n", | |
| "<br><br>\n", | |
| "ここで、変数間に完全な線形関係がある場合$$(X^TX)^{-1}$$の計算が不可能となる。<br>\n", | |
| "以下に3変数、3サンプルのデータを例に見ていく。<br>\n", | |
| "計画行列Xにおいて、変数$$x_1,x_2$$が次のように完全に線形関係であると想定する。($$x_2 = 3x_1$$)ここで、\n", | |
| "$$x_3$$は特に関連のない変数とする。\n", | |
| "<br><br>\n", | |
| "$$X = \\left(\n", | |
| " \\begin{array}{cccc}\n", | |
| " 1 & 3 & 3 \\\\\n", | |
| " 2 & 6 & 12 \\\\\n", | |
| " 5 & 15 & 5\n", | |
| " \\end{array}\n", | |
| " \\right)$$\n", | |
| "<br><br>\n", | |
| "$$X^T = \\left(\n", | |
| " \\begin{array}{cccc}\n", | |
| " 1 & 2 & 5 \\\\\n", | |
| " 3 & 6 & 15 \\\\\n", | |
| " 3 & 12 & 5\n", | |
| " \\end{array}\n", | |
| " \\right)$$\n", | |
| "<br><br>\n", | |
| "$$X^TX = \\left(\n", | |
| " \\begin{array}{cccc}\n", | |
| " 1 & 2 & 5 \\\\\n", | |
| " 3 & 6 & 15 \\\\\n", | |
| " 3 & 12 & 5\n", | |
| " \\end{array}\n", | |
| " \\right)\n", | |
| " \\left(\n", | |
| " \\begin{array}{cccc}\n", | |
| " 1 & 3 & 3 \\\\\n", | |
| " 2 & 6 & 12 \\\\\n", | |
| " 5 & 15 & 5\n", | |
| " \\end{array}\n", | |
| " \\right)\n", | |
| "= \\left(\n", | |
| " \\begin{array}{cccc}\n", | |
| " 30 & 90 & 52 \\\\\n", | |
| " 90 & 270 & 156 \\\\\n", | |
| " 52 & 156 & 178\n", | |
| " \\end{array}\n", | |
| " \\right) \n", | |
| "$$\n", | |
| "<br><br>\n", | |
| "ここで、3×3の行列Aの逆行列は以下のように求められる。\n", | |
| "<br><br>\n", | |
| "$$A=\\left(\n", | |
| " \\begin{array}{cccc}\n", | |
| " a_{11} & a_{12} & a_{13} \\\\\n", | |
| " a_{21} & a_{22} & a_{23} \\\\\n", | |
| " a_{31} & a_{32} & a_{33}\n", | |
| " \\end{array}\n", | |
| " \\right)$$\n", | |
| "<br>\n", | |
| "とおくと、\n", | |
| "<br>\n", | |
| "$$A^{-1}=\\frac{1}{detA} \\left(\n", | |
| " \\begin{array}{cccc}\n", | |
| " a_{22}a_{33}-a_{23}a_{32} & a_{13}a_{32}-a_{12}a_{33} & a_{12}a_{23}-a_{13}a_{22} \\\\\n", | |
| " a_{23}a_{31}-a_{21}a_{33} & a_{11}a_{33}-a_{13}a_{31} & a_{13}a_{21}-a_{11}a_{23} \\\\\n", | |
| " a_{21}a_{32}-a_{22}a_{31} & a_{12}a_{31}-a_{11}a_{32} & a_{11}a_{22}-a_{12}a_{21}\n", | |
| " \\end{array}\n", | |
| " \\right)$$\n", | |
| "<br>\n", | |
| "ここで、\n", | |
| "<br>\n", | |
| "$$detA = a_{11}a_{22}a_{33} + a_{21}a_{32}a_{13} + a_{31}a_{12}a_{23}\n", | |
| " - a_{11}a_{32}a_{23} - a_{31}a_{22}a_{13} - a_{21}a_{12}a_{33}$$\n", | |
| "<br><br>\n", | |
| "これを元に、先の計画行列のdetAを求めてみる。\n", | |
| "<br><br>\n", | |
| "$$detA = 30\\cdot 270\\cdot 178 + 90\\cdot 156\\cdot 52 + 52\\cdot 90\\cdot 156\n", | |
| " - 30\\cdot 156\\cdot 156 - 52\\cdot 270\\cdot 52 - 90\\cdot 90\\cdot 178$$\n", | |
| "$$= 1441800 + 730080 + 730080 - 730080 - 730080 - 1441800 = 0$$" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "このように、完全な多重共線性(変数間に完全な線形関係)があると逆行列計算の際にdetAが0(分母が0)となり計算ができなくなることがわかる。 \n", | |
| " \n", | |
| "ここで、**リッジ回帰**の場合を考えてみる。 \n", | |
| "リッジ回帰の計算については[ここ](https://gist.github.com/shotahorii/529a8c979ba0836a34c383c79c1e9f4b)参照。" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 67, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/latex": [ | |
| "リッジ回帰の解は以下のように表される。\n", | |
| "<br><br>\n", | |
| "$$\\beta_{ridge} = (X^TX+\\lambda I)^{-1}X^Ty$$<br>\n", | |
| "λは定数、Iは単位行列である。\n", | |
| "<br><br>\n", | |
| "ここで、今回の例で考える。(λ=3とする。)\n", | |
| "<br><br>\n", | |
| "$$X^TX = \\left(\n", | |
| " \\begin{array}{cccc}\n", | |
| " 30 & 90 & 52 \\\\\n", | |
| " 90 & 270 & 156 \\\\\n", | |
| " 52 & 156 & 178\n", | |
| " \\end{array}\n", | |
| " \\right) \n", | |
| "$$\n", | |
| "<br><br>\n", | |
| "$$\\lambda I = \\left(\n", | |
| " \\begin{array}{cccc}\n", | |
| " 3 & 0 & 0 \\\\\n", | |
| " 0 & 3 & 0 \\\\\n", | |
| " 0 & 0 & 3\n", | |
| " \\end{array}\n", | |
| " \\right) \n", | |
| "$$\n", | |
| "<br><br>\n", | |
| "$$X^TX + \\lambda I = \\left(\n", | |
| " \\begin{array}{cccc}\n", | |
| " 30 & 90 & 52 \\\\\n", | |
| " 90 & 270 & 156 \\\\\n", | |
| " 52 & 156 & 178\n", | |
| " \\end{array}\n", | |
| " \\right) + \n", | |
| " \\left(\n", | |
| " \\begin{array}{cccc}\n", | |
| " 3 & 0 & 0 \\\\\n", | |
| " 0 & 3 & 0 \\\\\n", | |
| " 0 & 0 & 3\n", | |
| " \\end{array}\n", | |
| " \\right) = \\left(\n", | |
| " \\begin{array}{cccc}\n", | |
| " 33 & 90 & 52 \\\\\n", | |
| " 90 & 273 & 156 \\\\\n", | |
| " 52 & 156 & 181\n", | |
| " \\end{array}\n", | |
| " \\right)\n", | |
| "$$\n", | |
| "<br><br>\n", | |
| "ここで、$$det(X^TX + \\lambda I)$$を計算する。\n", | |
| "<br><br>\n", | |
| "$$det(X^TX + \\lambda I) = 33\\cdot 273\\cdot 181 + 90\\cdot 156\\cdot 52 + 52\\cdot 90\\cdot 156\n", | |
| " - 33\\cdot 156\\cdot 156 - 52\\cdot 273\\cdot 52 - 90\\cdot 90\\cdot 181$$<br>\n", | |
| "$$= 1630629 + 730080 + 730080 - 803088 - 738192 - 1466100 = 83409$$\n", | |
| "<br><br>\n", | |
| "となり、多重共線性があってもdetAの分母が0にならず係数を計算できていることがわかる。<br>\n", | |
| "逆行列が存在する行列を「正則」という。<b>正則化</b>を行うことで多重共線性のために正則でなかった行列が正則となる。" | |
| ], | |
| "text/plain": [ | |
| "<IPython.core.display.Latex object>" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "%%latex\n", | |
| "リッジ回帰の解は以下のように表される。\n", | |
| "<br><br>\n", | |
| "$$\\beta_{ridge} = (X^TX+\\lambda I)^{-1}X^Ty$$<br>\n", | |
| "λは定数、Iは単位行列である。\n", | |
| "<br><br>\n", | |
| "ここで、今回の例で考える。(λ=3とする。)\n", | |
| "<br><br>\n", | |
| "$$X^TX = \\left(\n", | |
| " \\begin{array}{cccc}\n", | |
| " 30 & 90 & 52 \\\\\n", | |
| " 90 & 270 & 156 \\\\\n", | |
| " 52 & 156 & 178\n", | |
| " \\end{array}\n", | |
| " \\right) \n", | |
| "$$\n", | |
| "<br><br>\n", | |
| "$$\\lambda I = \\left(\n", | |
| " \\begin{array}{cccc}\n", | |
| " 3 & 0 & 0 \\\\\n", | |
| " 0 & 3 & 0 \\\\\n", | |
| " 0 & 0 & 3\n", | |
| " \\end{array}\n", | |
| " \\right) \n", | |
| "$$\n", | |
| "<br><br>\n", | |
| "$$X^TX + \\lambda I = \\left(\n", | |
| " \\begin{array}{cccc}\n", | |
| " 30 & 90 & 52 \\\\\n", | |
| " 90 & 270 & 156 \\\\\n", | |
| " 52 & 156 & 178\n", | |
| " \\end{array}\n", | |
| " \\right) + \n", | |
| " \\left(\n", | |
| " \\begin{array}{cccc}\n", | |
| " 3 & 0 & 0 \\\\\n", | |
| " 0 & 3 & 0 \\\\\n", | |
| " 0 & 0 & 3\n", | |
| " \\end{array}\n", | |
| " \\right) = \\left(\n", | |
| " \\begin{array}{cccc}\n", | |
| " 33 & 90 & 52 \\\\\n", | |
| " 90 & 273 & 156 \\\\\n", | |
| " 52 & 156 & 181\n", | |
| " \\end{array}\n", | |
| " \\right)\n", | |
| "$$\n", | |
| "<br><br>\n", | |
| "ここで、$$det(X^TX + \\lambda I)$$を計算する。\n", | |
| "<br><br>\n", | |
| "$$det(X^TX + \\lambda I) = 33\\cdot 273\\cdot 181 + 90\\cdot 156\\cdot 52 + 52\\cdot 90\\cdot 156\n", | |
| " - 33\\cdot 156\\cdot 156 - 52\\cdot 273\\cdot 52 - 90\\cdot 90\\cdot 181$$<br>\n", | |
| "$$= 1630629 + 730080 + 730080 - 803088 - 738192 - 1466100 = 83409$$\n", | |
| "<br><br>\n", | |
| "となり、多重共線性があってもdetAの分母が0にならず係数を計算できていることがわかる。<br>\n", | |
| "逆行列が存在する行列を「正則」という。<b>正則化</b>を行うことで多重共線性のために正則でなかった行列が正則となる。" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "次に、最初にあげた多重共線性の問題の2点目 - 変数間に完全な線形関係はないが強い相関があるケース - について考える。\n", | |
| "#### 2. 変数間に完全な線形関係はないが強い相関があるケース" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 71, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/latex": [ | |
| "上で見たように、通常の最小二乗法において変数間に完全な線形関係がある場合には$$det(X^TX)$$が0となり解が求まらない。\n", | |
| "一方、完全な線形関係はないが強い相関があるという場合$$det(X^TX)$$は0に近づき、逆行列が無限大へと発散する。\n", | |
| "すなわち、推定量の値が非常に不安定になる。\n", | |
| "<br><br>\n", | |
| "ここで、リッジ回帰は$$det(X^TX)$$の値を0から遠ざけることを考えると、この完全な線形関係はないが強い相関がある\n", | |
| "ケースについても有効であることがわかる。\n", | |
| "<br><br>\n", | |
| "ただし、説明変数が応答変数に対してどの程度影響を与えているか、という推定量の信頼性について、どの程度のλを設定\n", | |
| "すれば信頼できると言えるか、というのは実際わからない。<br>\n", | |
| "ただこれについては、実務において、説明変数の応答変数に対する影響度を信頼できる値として定量的に示す必要があるような\n", | |
| "特殊な場合を除き、Cross Validationなどで「当てはまりのよいモデルかつ予測能力も高い」ことが示されていれば\n", | |
| "問題ないように思う。" | |
| ], | |
| "text/plain": [ | |
| "<IPython.core.display.Latex object>" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "%%latex\n", | |
| "上で見たように、通常の最小二乗法において変数間に完全な線形関係がある場合には$$det(X^TX)$$が0となり解が求まらない。\n", | |
| "一方、完全な線形関係はないが強い相関があるという場合$$det(X^TX)$$は0に近づき、逆行列が無限大へと発散する。\n", | |
| "すなわち、推定量の値が非常に不安定になる。\n", | |
| "<br><br>\n", | |
| "ここで、リッジ回帰は$$det(X^TX)$$の値を0から遠ざけることを考えると、この完全な線形関係はないが強い相関がある\n", | |
| "ケースについても有効であることがわかる。\n", | |
| "<br><br>\n", | |
| "ただし、説明変数が応答変数に対してどの程度影響を与えているか、という推定量の信頼性について、どの程度のλを設定\n", | |
| "すれば信頼できると言えるか、というのは実際わからない。<br>\n", | |
| "ただこれについては、実務において、説明変数の応答変数に対する影響度を信頼できる値として定量的に示す必要があるような\n", | |
| "特殊な場合を除き、Cross Validationなどで「当てはまりのよいモデルかつ予測能力も高い」ことが示されていれば\n", | |
| "問題ないように思う。" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "#### データで確認\n", | |
| "完全ではないが強い相関を持つ変数を持つデータを用いて、学習に使うデータサンプルを少しずつ変え、推定量がどの程度安定しているかOLSとRidge回帰を比較してみる。" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 145, | |
| "metadata": { | |
| "collapsed": true | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "import math\n", | |
| "import random\n", | |
| "import numpy as np\n", | |
| "import pandas as pd\n", | |
| "from functools import reduce\n", | |
| "\n", | |
| "from sklearn.datasets import load_boston\n", | |
| "\n", | |
| "from matplotlib import pyplot as plt\n", | |
| "%matplotlib inline" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "##### Models" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 273, | |
| "metadata": { | |
| "collapsed": true | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "def fitting_ols(X,y):\n", | |
| " beta = np.dot(np.dot(np.linalg.inv(np.dot(X.T,X)),X.T),y) # raw implementation\n", | |
| " #(beta, residuals, rank, s) = np.linalg.lstsq(X, y) # numpy function for least square \n", | |
| " return beta\n", | |
| "\n", | |
| "def fitting_ridge(X,y,l):\n", | |
| " beta = np.dot(np.dot(np.linalg.inv(np.dot(X.T,X)+l*np.eye(len(X.T))),X.T),y) # raw implementation\n", | |
| " return beta" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "##### Data" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 274, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "d = load_boston()\n", | |
| "y = d.target\n", | |
| "X = d.data" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 275, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/html": [ | |
| "<div>\n", | |
| "<table border=\"1\" class=\"dataframe\">\n", | |
| " <thead>\n", | |
| " <tr style=\"text-align: right;\">\n", | |
| " <th></th>\n", | |
| " <th>x1</th>\n", | |
| " <th>x2</th>\n", | |
| " <th>x3</th>\n", | |
| " <th>x4</th>\n", | |
| " <th>x5</th>\n", | |
| " <th>x6</th>\n", | |
| " <th>x7</th>\n", | |
| " <th>x8</th>\n", | |
| " <th>x9</th>\n", | |
| " <th>x10</th>\n", | |
| " <th>x11</th>\n", | |
| " <th>x12</th>\n", | |
| " <th>x13</th>\n", | |
| " <th>y</th>\n", | |
| " </tr>\n", | |
| " </thead>\n", | |
| " <tbody>\n", | |
| " <tr>\n", | |
| " <th>x1</th>\n", | |
| " <td>1.000000</td>\n", | |
| " <td>-0.200469</td>\n", | |
| " <td>0.406583</td>\n", | |
| " <td>-0.055892</td>\n", | |
| " <td>0.420972</td>\n", | |
| " <td>-0.219247</td>\n", | |
| " <td>0.352734</td>\n", | |
| " <td>-0.379670</td>\n", | |
| " <td>0.625505</td>\n", | |
| " <td>0.582764</td>\n", | |
| " <td>0.289946</td>\n", | |
| " <td>-0.385064</td>\n", | |
| " <td>0.455621</td>\n", | |
| " <td>-0.388305</td>\n", | |
| " </tr>\n", | |
| " <tr>\n", | |
| " <th>x2</th>\n", | |
| " <td>-0.200469</td>\n", | |
| " <td>1.000000</td>\n", | |
| " <td>-0.533828</td>\n", | |
| " <td>-0.042697</td>\n", | |
| " <td>-0.516604</td>\n", | |
| " <td>0.311991</td>\n", | |
| " <td>-0.569537</td>\n", | |
| " <td>0.664408</td>\n", | |
| " <td>-0.311948</td>\n", | |
| " <td>-0.314563</td>\n", | |
| " <td>-0.391679</td>\n", | |
| " <td>0.175520</td>\n", | |
| " <td>-0.412995</td>\n", | |
| " <td>0.360445</td>\n", | |
| " </tr>\n", | |
| " <tr>\n", | |
| " <th>x3</th>\n", | |
| " <td>0.406583</td>\n", | |
| " <td>-0.533828</td>\n", | |
| " <td>1.000000</td>\n", | |
| " <td>0.062938</td>\n", | |
| " <td>0.763651</td>\n", | |
| " <td>-0.391676</td>\n", | |
| " <td>0.644779</td>\n", | |
| " <td>-0.708027</td>\n", | |
| " <td>0.595129</td>\n", | |
| " <td>0.720760</td>\n", | |
| " <td>0.383248</td>\n", | |
| " <td>-0.356977</td>\n", | |
| " <td>0.603800</td>\n", | |
| " <td>-0.483725</td>\n", | |
| " </tr>\n", | |
| " <tr>\n", | |
| " <th>x4</th>\n", | |
| " <td>-0.055892</td>\n", | |
| " <td>-0.042697</td>\n", | |
| " <td>0.062938</td>\n", | |
| " <td>1.000000</td>\n", | |
| " <td>0.091203</td>\n", | |
| " <td>0.091251</td>\n", | |
| " <td>0.086518</td>\n", | |
| " <td>-0.099176</td>\n", | |
| " <td>-0.007368</td>\n", | |
| " <td>-0.035587</td>\n", | |
| " <td>-0.121515</td>\n", | |
| " <td>0.048788</td>\n", | |
| " <td>-0.053929</td>\n", | |
| " <td>0.175260</td>\n", | |
| " </tr>\n", | |
| " <tr>\n", | |
| " <th>x5</th>\n", | |
| " <td>0.420972</td>\n", | |
| " <td>-0.516604</td>\n", | |
| " <td>0.763651</td>\n", | |
| " <td>0.091203</td>\n", | |
| " <td>1.000000</td>\n", | |
| " <td>-0.302188</td>\n", | |
| " <td>0.731470</td>\n", | |
| " <td>-0.769230</td>\n", | |
| " <td>0.611441</td>\n", | |
| " <td>0.668023</td>\n", | |
| " <td>0.188933</td>\n", | |
| " <td>-0.380051</td>\n", | |
| " <td>0.590879</td>\n", | |
| " <td>-0.427321</td>\n", | |
| " </tr>\n", | |
| " <tr>\n", | |
| " <th>x6</th>\n", | |
| " <td>-0.219247</td>\n", | |
| " <td>0.311991</td>\n", | |
| " <td>-0.391676</td>\n", | |
| " <td>0.091251</td>\n", | |
| " <td>-0.302188</td>\n", | |
| " <td>1.000000</td>\n", | |
| " <td>-0.240265</td>\n", | |
| " <td>0.205246</td>\n", | |
| " <td>-0.209847</td>\n", | |
| " <td>-0.292048</td>\n", | |
| " <td>-0.355501</td>\n", | |
| " <td>0.128069</td>\n", | |
| " <td>-0.613808</td>\n", | |
| " <td>0.695360</td>\n", | |
| " </tr>\n", | |
| " <tr>\n", | |
| " <th>x7</th>\n", | |
| " <td>0.352734</td>\n", | |
| " <td>-0.569537</td>\n", | |
| " <td>0.644779</td>\n", | |
| " <td>0.086518</td>\n", | |
| " <td>0.731470</td>\n", | |
| " <td>-0.240265</td>\n", | |
| " <td>1.000000</td>\n", | |
| " <td>-0.747881</td>\n", | |
| " <td>0.456022</td>\n", | |
| " <td>0.506456</td>\n", | |
| " <td>0.261515</td>\n", | |
| " <td>-0.273534</td>\n", | |
| " <td>0.602339</td>\n", | |
| " <td>-0.376955</td>\n", | |
| " </tr>\n", | |
| " <tr>\n", | |
| " <th>x8</th>\n", | |
| " <td>-0.379670</td>\n", | |
| " <td>0.664408</td>\n", | |
| " <td>-0.708027</td>\n", | |
| " <td>-0.099176</td>\n", | |
| " <td>-0.769230</td>\n", | |
| " <td>0.205246</td>\n", | |
| " <td>-0.747881</td>\n", | |
| " <td>1.000000</td>\n", | |
| " <td>-0.494588</td>\n", | |
| " <td>-0.534432</td>\n", | |
| " <td>-0.232471</td>\n", | |
| " <td>0.291512</td>\n", | |
| " <td>-0.496996</td>\n", | |
| " <td>0.249929</td>\n", | |
| " </tr>\n", | |
| " <tr>\n", | |
| " <th>x9</th>\n", | |
| " <td>0.625505</td>\n", | |
| " <td>-0.311948</td>\n", | |
| " <td>0.595129</td>\n", | |
| " <td>-0.007368</td>\n", | |
| " <td>0.611441</td>\n", | |
| " <td>-0.209847</td>\n", | |
| " <td>0.456022</td>\n", | |
| " <td>-0.494588</td>\n", | |
| " <td>1.000000</td>\n", | |
| " <td>0.910228</td>\n", | |
| " <td>0.464741</td>\n", | |
| " <td>-0.444413</td>\n", | |
| " <td>0.488676</td>\n", | |
| " <td>-0.381626</td>\n", | |
| " </tr>\n", | |
| " <tr>\n", | |
| " <th>x10</th>\n", | |
| " <td>0.582764</td>\n", | |
| " <td>-0.314563</td>\n", | |
| " <td>0.720760</td>\n", | |
| " <td>-0.035587</td>\n", | |
| " <td>0.668023</td>\n", | |
| " <td>-0.292048</td>\n", | |
| " <td>0.506456</td>\n", | |
| " <td>-0.534432</td>\n", | |
| " <td>0.910228</td>\n", | |
| " <td>1.000000</td>\n", | |
| " <td>0.460853</td>\n", | |
| " <td>-0.441808</td>\n", | |
| " <td>0.543993</td>\n", | |
| " <td>-0.468536</td>\n", | |
| " </tr>\n", | |
| " <tr>\n", | |
| " <th>x11</th>\n", | |
| " <td>0.289946</td>\n", | |
| " <td>-0.391679</td>\n", | |
| " <td>0.383248</td>\n", | |
| " <td>-0.121515</td>\n", | |
| " <td>0.188933</td>\n", | |
| " <td>-0.355501</td>\n", | |
| " <td>0.261515</td>\n", | |
| " <td>-0.232471</td>\n", | |
| " <td>0.464741</td>\n", | |
| " <td>0.460853</td>\n", | |
| " <td>1.000000</td>\n", | |
| " <td>-0.177383</td>\n", | |
| " <td>0.374044</td>\n", | |
| " <td>-0.507787</td>\n", | |
| " </tr>\n", | |
| " <tr>\n", | |
| " <th>x12</th>\n", | |
| " <td>-0.385064</td>\n", | |
| " <td>0.175520</td>\n", | |
| " <td>-0.356977</td>\n", | |
| " <td>0.048788</td>\n", | |
| " <td>-0.380051</td>\n", | |
| " <td>0.128069</td>\n", | |
| " <td>-0.273534</td>\n", | |
| " <td>0.291512</td>\n", | |
| " <td>-0.444413</td>\n", | |
| " <td>-0.441808</td>\n", | |
| " <td>-0.177383</td>\n", | |
| " <td>1.000000</td>\n", | |
| " <td>-0.366087</td>\n", | |
| " <td>0.333461</td>\n", | |
| " </tr>\n", | |
| " <tr>\n", | |
| " <th>x13</th>\n", | |
| " <td>0.455621</td>\n", | |
| " <td>-0.412995</td>\n", | |
| " <td>0.603800</td>\n", | |
| " <td>-0.053929</td>\n", | |
| " <td>0.590879</td>\n", | |
| " <td>-0.613808</td>\n", | |
| " <td>0.602339</td>\n", | |
| " <td>-0.496996</td>\n", | |
| " <td>0.488676</td>\n", | |
| " <td>0.543993</td>\n", | |
| " <td>0.374044</td>\n", | |
| " <td>-0.366087</td>\n", | |
| " <td>1.000000</td>\n", | |
| " <td>-0.737663</td>\n", | |
| " </tr>\n", | |
| " <tr>\n", | |
| " <th>y</th>\n", | |
| " <td>-0.388305</td>\n", | |
| " <td>0.360445</td>\n", | |
| " <td>-0.483725</td>\n", | |
| " <td>0.175260</td>\n", | |
| " <td>-0.427321</td>\n", | |
| " <td>0.695360</td>\n", | |
| " <td>-0.376955</td>\n", | |
| " <td>0.249929</td>\n", | |
| " <td>-0.381626</td>\n", | |
| " <td>-0.468536</td>\n", | |
| " <td>-0.507787</td>\n", | |
| " <td>0.333461</td>\n", | |
| " <td>-0.737663</td>\n", | |
| " <td>1.000000</td>\n", | |
| " </tr>\n", | |
| " </tbody>\n", | |
| "</table>\n", | |
| "</div>" | |
| ], | |
| "text/plain": [ | |
| " x1 x2 x3 x4 x5 x6 x7 \\\n", | |
| "x1 1.000000 -0.200469 0.406583 -0.055892 0.420972 -0.219247 0.352734 \n", | |
| "x2 -0.200469 1.000000 -0.533828 -0.042697 -0.516604 0.311991 -0.569537 \n", | |
| "x3 0.406583 -0.533828 1.000000 0.062938 0.763651 -0.391676 0.644779 \n", | |
| "x4 -0.055892 -0.042697 0.062938 1.000000 0.091203 0.091251 0.086518 \n", | |
| "x5 0.420972 -0.516604 0.763651 0.091203 1.000000 -0.302188 0.731470 \n", | |
| "x6 -0.219247 0.311991 -0.391676 0.091251 -0.302188 1.000000 -0.240265 \n", | |
| "x7 0.352734 -0.569537 0.644779 0.086518 0.731470 -0.240265 1.000000 \n", | |
| "x8 -0.379670 0.664408 -0.708027 -0.099176 -0.769230 0.205246 -0.747881 \n", | |
| "x9 0.625505 -0.311948 0.595129 -0.007368 0.611441 -0.209847 0.456022 \n", | |
| "x10 0.582764 -0.314563 0.720760 -0.035587 0.668023 -0.292048 0.506456 \n", | |
| "x11 0.289946 -0.391679 0.383248 -0.121515 0.188933 -0.355501 0.261515 \n", | |
| "x12 -0.385064 0.175520 -0.356977 0.048788 -0.380051 0.128069 -0.273534 \n", | |
| "x13 0.455621 -0.412995 0.603800 -0.053929 0.590879 -0.613808 0.602339 \n", | |
| "y -0.388305 0.360445 -0.483725 0.175260 -0.427321 0.695360 -0.376955 \n", | |
| "\n", | |
| " x8 x9 x10 x11 x12 x13 y \n", | |
| "x1 -0.379670 0.625505 0.582764 0.289946 -0.385064 0.455621 -0.388305 \n", | |
| "x2 0.664408 -0.311948 -0.314563 -0.391679 0.175520 -0.412995 0.360445 \n", | |
| "x3 -0.708027 0.595129 0.720760 0.383248 -0.356977 0.603800 -0.483725 \n", | |
| "x4 -0.099176 -0.007368 -0.035587 -0.121515 0.048788 -0.053929 0.175260 \n", | |
| "x5 -0.769230 0.611441 0.668023 0.188933 -0.380051 0.590879 -0.427321 \n", | |
| "x6 0.205246 -0.209847 -0.292048 -0.355501 0.128069 -0.613808 0.695360 \n", | |
| "x7 -0.747881 0.456022 0.506456 0.261515 -0.273534 0.602339 -0.376955 \n", | |
| "x8 1.000000 -0.494588 -0.534432 -0.232471 0.291512 -0.496996 0.249929 \n", | |
| "x9 -0.494588 1.000000 0.910228 0.464741 -0.444413 0.488676 -0.381626 \n", | |
| "x10 -0.534432 0.910228 1.000000 0.460853 -0.441808 0.543993 -0.468536 \n", | |
| "x11 -0.232471 0.464741 0.460853 1.000000 -0.177383 0.374044 -0.507787 \n", | |
| "x12 0.291512 -0.444413 -0.441808 -0.177383 1.000000 -0.366087 0.333461 \n", | |
| "x13 -0.496996 0.488676 0.543993 0.374044 -0.366087 1.000000 -0.737663 \n", | |
| "y 0.249929 -0.381626 -0.468536 -0.507787 0.333461 -0.737663 1.000000 " | |
| ] | |
| }, | |
| "execution_count": 275, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| } | |
| ], | |
| "source": [ | |
| "# check correlations between x variables and y\n", | |
| "cols = ['x'+str(i) for i in range(1,len(X.T)+1)]\n", | |
| "df = pd.DataFrame(X,columns=cols)\n", | |
| "df['y'] = y\n", | |
| "df.corr()" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 289, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "# use x6 and X8, then make a variable highly correlating with x6\n", | |
| "# original cols (x6 and x8)\n", | |
| "cond = [5,7]\n", | |
| "orig_X = X[:,cond]\n", | |
| "\n", | |
| "# making correlated variables\n", | |
| "corr_x1 = X[:,5] + 0.1*np.random.randn(len(X))\n", | |
| "corr_x2 = X[:,7] + 0.1*np.random.randn(len(X))\n", | |
| "\n", | |
| "# merge\n", | |
| "multico_X = np.insert(orig_X, 2, corr_x1, axis=1)\n", | |
| "multico_X = np.insert(multico_X, 3, corr_x2, axis=1)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 277, | |
| "metadata": { | |
| "collapsed": true | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "# dataset creation functions\n", | |
| "def add_intercept(X):\n", | |
| " return np.insert(X, 0, np.ones(len(X)), axis=1)\n", | |
| "\n", | |
| "def dataset(X,y,p):\n", | |
| " n = round(len(X)*p)\n", | |
| " ind = sorted(random.sample(range(0,len(X)),n))\n", | |
| " return add_intercept(X[ind,:]), y[ind]" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "##### Main" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 308, | |
| "metadata": { | |
| "collapsed": true | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "p = 0.95 # use 90% of original samples\n", | |
| "itr = 1000 # number of trial" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 318, | |
| "metadata": { | |
| "collapsed": true | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "# No multicollenearity\n", | |
| "ols_betas = []\n", | |
| "ridge_betas = []\n", | |
| "for i in range(itr):\n", | |
| " sub_X, sub_y = dataset(orig_X,y,p)\n", | |
| " ols_betas.append(fitting_ols(sub_X,sub_y))\n", | |
| " ridge_betas.append(fitting_ridge(sub_X,sub_y,3)) # lambda = 3" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 319, | |
| "metadata": { | |
| "collapsed": true | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "df_ols_orig = pd.DataFrame(ols_betas,columns=['b0','b1','b2'])\n", | |
| "df_ridge_orig = pd.DataFrame(ridge_betas,columns=['b0','b1','b2'])" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 320, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": [ | |
| "b0 1.011566\n", | |
| "b1 0.026488\n", | |
| "b2 0.001037\n", | |
| "dtype: float64" | |
| ] | |
| }, | |
| "execution_count": 320, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| } | |
| ], | |
| "source": [ | |
| "df_ols_orig.var()" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 321, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": [ | |
| "b0 0.410973\n", | |
| "b1 0.012125\n", | |
| "b2 0.001085\n", | |
| "dtype: float64" | |
| ] | |
| }, | |
| "execution_count": 321, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| } | |
| ], | |
| "source": [ | |
| "df_ridge_orig.var()" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 322, | |
| "metadata": { | |
| "collapsed": true | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "ols_betas = []\n", | |
| "ridge_betas = []\n", | |
| "for i in range(itr):\n", | |
| " sub_X, sub_y = dataset(multico_X,y,p)\n", | |
| " ols_betas.append(fitting_ols(sub_X,sub_y))\n", | |
| " ridge_betas.append(fitting_ridge(sub_X,sub_y,3)) # lambda = 3" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 323, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "df_ols_multico = pd.DataFrame(ols_betas,columns=['b0','b1','b2','b3','b4'])\n", | |
| "df_ridge_multico = pd.DataFrame(ridge_betas,columns=['b0','b1','b2','b3','b4'])" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 324, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": [ | |
| "b0 1.014518\n", | |
| "b1 0.428334\n", | |
| "b2 0.479558\n", | |
| "b3 0.406282\n", | |
| "b4 0.482585\n", | |
| "dtype: float64" | |
| ] | |
| }, | |
| "execution_count": 324, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| } | |
| ], | |
| "source": [ | |
| "df_ols_multico.var()" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 325, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": [ | |
| "b0 0.407626\n", | |
| "b1 0.082622\n", | |
| "b2 0.100647\n", | |
| "b3 0.085184\n", | |
| "b4 0.102217\n", | |
| "dtype: float64" | |
| ] | |
| }, | |
| "execution_count": 325, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| } | |
| ], | |
| "source": [ | |
| "df_ridge_multico.var()" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "OLSでは、多重共線性があるデータを用いた場合に推定量の分散が大きくなっていることが見て取れる。 \n", | |
| "一方で、λ=3としたリッジ回帰では推定量の分散は小さく抑えられている。ただし、λを大きく取ればバイアスは大きくなる点に注意。" | |
| ] | |
| } | |
| ], | |
| "metadata": { | |
| "kernelspec": { | |
| "display_name": "Python 3", | |
| "language": "python", | |
| "name": "python3" | |
| }, | |
| "language_info": { | |
| "codemirror_mode": { | |
| "name": "ipython", | |
| "version": 3 | |
| }, | |
| "file_extension": ".py", | |
| "mimetype": "text/x-python", | |
| "name": "python", | |
| "nbconvert_exporter": "python", | |
| "pygments_lexer": "ipython3", | |
| "version": "3.5.2" | |
| } | |
| }, | |
| "nbformat": 4, | |
| "nbformat_minor": 0 | |
| } |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment