ex09notebookA.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "view-in-github"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/takatakamanbou/Data/blob/2021/ex09notebookA.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2UzDNLmVS3po"
      },
      "source": [
        "# Data2021 ex09notebookA\n",
        "\n",
        "<img width=64 src=\"https://www-tlab.math.ryukoku.ac.jp/~takataka/course/Data/Data-logo09.png\"> https://www-tlab.math.ryukoku.ac.jp/wiki/?Data/2021\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7H45eNIcd2xB"
      },
      "source": [
        "----\n",
        "## 回帰分析入門(3) \n",
        "----"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WM4RlkyhP5yo"
      },
      "source": [
        "----\n",
        "### 最小二乗法による直線あてはめの解の導出\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "hQd6yDT-di3D"
      },
      "source": [
        "2次元のデータに直線を当てはめる最小二乗法の問題は，次のように表されます．\n",
        "\n",
        "> データ  $(x_n, y_n)$ （$n = 1, 2, \\dots , N$） が与えられたときに，\n",
        ">\n",
        "> $$\n",
        "E(\\hat{a}, \\hat{b}) = \\sum_{n=1}^{N}(y_n - (\\hat{a}x_n + \\hat{b}))^2\n",
        "$$\n",
        ">\n",
        "> を最小にする $(\\hat{a}, \\hat{b})$ を求めよ．\n",
        "\n",
        "すでに，この問題の解が次のようになることは説明しました．\n",
        "\n",
        "> $ E(\\hat{a}, \\hat{b}) $ を最小にする $(\\hat{a}, \\hat{b})$ は，次式を満たす．\n",
        ">\n",
        "> $$\n",
        "\\hat{a} = \\frac{s_{xy}}{s_x^2}\\qquad \\hat{b} = \\bar{y} - \\hat{a}\\bar{x} \\tag{*}\n",
        "$$\n",
        "> \n",
        "> ただし， $\\bar{x}, \\bar{y}$ は $x, y$ の平均，$s_{x}^2$ は $x$ の分散，$s_{xy}$ は $x$ と $y$ の共分散，つまり，\n",
        "> \n",
        "> $$\n",
        "\\begin{aligned}\n",
        "\\bar{x} &= \\frac{1}{N}\\sum_{n=1}^N x_n\\\\\n",
        "\\bar{y} &= \\frac{1}{N}\\sum_{n=1}^N y_n\\\\\n",
        "s_{x}^2 &= \\frac{1}{N}\\sum_{n=1}^{N}(x_n - \\bar{x})^2\\\\\n",
        "s_{xy} &= \\frac{1}{N}\\sum_{n=1}^{N}(x_n - \\bar{x})(y_n - \\bar{y})\\\\\n",
        "\\end{aligned}\n",
        "$$\n",
        "\n",
        "ここでは，このことを証明してみましょう．\n",
        "\n",
        "以下，簡単のため，回帰係数の推定値の $\\hat{}$ を省略して $a, b$ と表記していることがあります．\n",
        "また，$\\Sigma$の上下の式を省略していることがあります．"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "H0r86M-fTW8T"
      },
      "source": [
        "方針: $E(a, b)$ は， 2つの変数 $a, b$ に対する2次式であり，いわゆる「下に凸」な形なので，「$E(a, b)$ を平方完成する」か，「$E(a, b)$ を $a, b$ で偏微分して $0$ とおく」ことで，最小解を得ることができます．\n",
        "ここでは後者のやり方で考えてみます．"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zQD4QuJ-7Q9Q"
      },
      "source": [
        "$E(a, b)$ を $a$ で偏微分すると，\n",
        "\n",
        "$$\n",
        "\\begin{aligned}\n",
        "\\frac{\\partial E(a, b)}{\\partial a} &= \\frac{\\partial}{\\partial a} \\sum_{n=1}^{N}(y_n - (ax_n + b))^2 \\\\\n",
        "&= \\sum_{n} \\frac{\\partial}{\\partial a}(y_n - (ax_n + b))^2 \\\\\n",
        "&= \\sum_{n} 2(y_n - (ax_n + b)) \\frac{\\partial}{\\partial a}(y_n - (ax_n + b)) \\\\\n",
        "&= \\sum_{n} 2(y_n - (ax_n + b))(-x_n)\\\\\n",
        "&= 2\\sum_{n}(ax_n^2+bx_n - x_ny_n)\n",
        "\\end{aligned}\n",
        "$$\n",
        "\n",
        "となります．\n",
        "したがって，$\\frac{\\partial E(a, b)}{\\partial a} = 0$ とおくと，次の方程式が得られます．\n",
        "\n",
        "$$\n",
        "\\sum_{n=1}^{N}(ax_n^2+bx_n - x_ny_n) = a\\sum_{n=1}^{N} x_n^2 + b\\sum_{n=1}^{N} x_n - \\sum_{n=1}^{N}x_ny_n = 0 \\tag{1}\n",
        "$$\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xvubfbeJ9g01"
      },
      "source": [
        "\n",
        "また， $E(a, b)$ を $b$ で偏微分すると，\n",
        "\n",
        "$$\n",
        "\\begin{aligned}\n",
        "\\frac{\\partial E(a, b)}{\\partial b} &= \\frac{\\partial}{\\partial b} \\sum_{n=1}^{N}(y_n - (ax_n + b))^2 \\\\\n",
        "&= \\sum_{n} \\frac{\\partial}{\\partial b}(y_n - (ax_n + b))^2 \\\\\n",
        "&= \\sum_{n} 2(y_n - (ax_n + b)) \\frac{\\partial}{\\partial b}(y_n - (ax_n + b)) \\\\\n",
        "&= \\sum_{n} 2(y_n - (ax_n + b))(-1)\\\\\n",
        "&= 2\\sum_{n=1}^{N}(ax_n + b - y_n)\n",
        "\\end{aligned}\n",
        "$$\n",
        "\n",
        "となります．したがって，$\\frac{\\partial E(a, b)}{\\partial b} = 0$ とおくと，次の方程式が得られます．\n",
        "\n",
        "$$\n",
        "\\sum_{n=1}^{N}(ax_n + b - y_n) =  a\\sum_{n=1}^{N}x_n + b\\sum_{n=1}^{N}1 - \\sum_{n=1}^{N}y_n = 0 \\tag{2}\n",
        "$$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "FBPsIXCtXCEv"
      },
      "source": [
        "以上より，$E(a, b)$ を最小にする $(a, b)$ は，式(1)と(2)の二つの方程式を連立した方程式の解でなければならないことがわかります．\n",
        "それぞれの式の左辺の最後の項を右辺へ移項した形で書くと，この連立方程式は次のようになります．\n",
        "\n",
        "$$\n",
        "\\begin{aligned}\n",
        "a\\sum_{n=1}^{N}x_n^2 + b\\sum_{n=1}^{N}x_n &= \\sum_{n=1}^{N}x_ny_n \\qquad (3)\\\\\n",
        "a\\sum_{n=1}^{N}x_n + b\\sum_{n=1}^{N}1 &= \\sum_{n=1}^{N}y_n\\quad\\qquad (4)\n",
        "\\end{aligned}\n",
        "$$\n",
        "\n",
        "この連立方程式のことを，**正規方程式** と呼ぶことがあります．参考までに，行列を使って表したものも示しておきます．\n",
        "\n",
        "$$\n",
        "\\begin{pmatrix}\n",
        "\\displaystyle\\sum_{n=1}^{N}x_n^2 & \\displaystyle\\sum_{n=1}^{N}x_n\\\\\n",
        "\\displaystyle\\sum_{n=1}^{N}x_n & \\displaystyle\\sum_{n=1}^{N}1\\\\\n",
        "\\end{pmatrix}\n",
        "\\begin{pmatrix}\n",
        "a \\\\ b\n",
        "\\end{pmatrix} =\n",
        "\\begin{pmatrix}\n",
        "\\displaystyle\\sum_{n=1}^{N}x_ny_n\\\\\n",
        "\\displaystyle\\sum_{n=1}^{N}y_n\n",
        "\\end{pmatrix} \n",
        "$$\n",
        "\n",
        "\n",
        "正規方程式を $a, b$ について解けば，$E(a, b)$ を最小にする $(a, b)$ が求まります．\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "clH1NocJQ6nl"
      },
      "source": [
        "まず，式(4)を $b$ について解いてみます．\n",
        "\n",
        "$$\n",
        "b\\times N = \\sum_{n=1}^Ny_n - a\\sum_{n=1}^{N}x_n\n",
        "$$\n",
        "\n",
        "両辺を $N$ で割ると\n",
        "\n",
        "$$\n",
        "b = \\frac{1}{N}\\sum_{n=1}^{N}y_n - a\\times\\frac{1}{N}\\sum_{n=1}^{N}x_n = \\bar{y} - a\\bar{x}\n",
        "$$\n",
        "\n",
        "が得られます．これは，すでに示した正規方程式の解のうち，$\\hat{b}$の方の式に一致していますね（式(*)参照）．\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "hzQwT790N4NV"
      },
      "source": [
        "次に，この式を式(3)の$b$に代入して変形すると，\n",
        "\n",
        "$$\n",
        "\\begin{aligned}\n",
        "a\\sum_{n=1}^{N}x_n^2 + (\\bar{y} -a\\bar{x})\\sum_{n=1}^{N}x_n &= \n",
        "\\sum_{n=1}^{N}x_ny_n \\\\\n",
        "a\\sum_{n=1}^{N}x_n^2 - a\\bar{x}\\sum_{n=1}^{N}x_n+ \\bar{y} \\sum_{n=1}^{N}x_n &= \n",
        "\\sum_{n=1}^{N}x_ny_n \\\\\n",
        "a\\sum_{n=1}^{N}x_n^2 - a\\bar{x}\\sum_{n=1}^{N}x_n &= \n",
        "\\sum_{n=1}^{N}x_ny_n - \\bar{y} \\sum_{n=1}^{N}x_n\\\\\n",
        "\\end{aligned}\n",
        "$$\n",
        "\n",
        "となります．最後の等式の両辺を $N$ で割ると，\n",
        "\n",
        "$$\n",
        "\\begin{aligned}\n",
        "a\\times\\frac{1}{N}\\sum_{n=1}^{N}x_n^2 - a\\bar{x}\\times \\frac{1}{N}\\sum_{n=1}^{N}x_n &= \\frac{1}{N}\n",
        "\\sum_{n=1}^{N}x_ny_n - \\bar{y} \\times \\frac{1}{N}\\sum_{n=1}^{N}x_n\\\\\n",
        "a\\times\\left(\\frac{1}{N}\\sum_{n=1}^{N}x_n^2 - \\bar{x}^2 \\right) &= \\frac{1}{N}\n",
        "\\sum_{n=1}^{N}x_ny_n - \\bar{x}\\bar{y} \\qquad (5)\n",
        "\\end{aligned}\n",
        "$$\n",
        "\n",
        "となります．"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "oJDAP1LwQxoR"
      },
      "source": [
        "ここで，\n",
        "\n",
        "$$\n",
        "\\begin{aligned}\n",
        "s_x^2 &= \\frac{1}{N}\\sum_{n=1}^{N}(x_n - \\bar{x})^2 = \\frac{1}{N}\\sum_{n=1}^{N}(x_n^2 - 2\\bar{x}x_n + \\bar{x}^2)\\\\\n",
        "&= \\frac{1}{N}\\sum_{n=1}^{N}x_n^2 - 2\\bar{x}\\times\\frac{1}{N}\\sum_{n=1}^{N}x_n + \\bar{x}^2\\times\\frac{1}{N}\\sum_{n=1}^{N}1\\\\\n",
        "&= \\frac{1}{N}\\sum_{n=1}^{N}x_n^2 - 2\\bar{x}^2 + \\bar{x}^2 = \\frac{1}{N}\\sum_{n=1}^{N}x_n^2 - \\bar{x}^2\n",
        "\\end{aligned}\n",
        "$$\n",
        "\n",
        "および\n",
        "\n",
        "$$\n",
        "\\begin{aligned}\n",
        "s_{xy} &= \\frac{1}{N}\\sum_{n=1}^{N}(x_n - \\bar{x})(y_n-\\bar{y}) = \\frac{1}{N}\\sum_{n=1}^{N}(x_ny_n - \\bar{x}y_n - x_n\\bar{y} + \\bar{x}\\bar{y})\\\\\n",
        "&= \\frac{1}{N}\\sum_{n=1}^{N}x_ny_n - \\bar{x}\\times\\frac{1}{N}\\sum_{n=1}^{N}y_n - \\bar{y}\\times\\frac{1}{N}\\sum_{n=1}^{N}x_n+ \\bar{x}\\bar{y}\\times\\frac{1}{N}\\sum_{n=1}^{N}1\\\\\n",
        "&= \\frac{1}{N}\\sum_{n=1}^{N}x_ny_n - \\bar{x}\\bar{y}- \\bar{x}\\bar{y} + \\bar{x}\\bar{y} = \\frac{1}{N}\\sum_{n=1}^{N}x_ny_n - \\bar{x}\\bar{y}\n",
        "\\end{aligned}\n",
        "$$\n",
        "\n",
        "であることから，式(5)は次のように書き直せます．\n",
        "\n",
        "$$\n",
        "as_{x}^2 = s_{xy}\n",
        "$$\n",
        "\n",
        "これを $a$ について解けば，解 $a$ の式\n",
        "\n",
        "$$\n",
        "a = \\frac{s_{xy}}{s_x^2}\n",
        "$$\n",
        "\n",
        "が得られます．\n",
        "\n",
        "以上より，$ E(a, b) $ を最小にする $(a, b)$ は，式(*)で求まることが示されました．\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dogmFuo-dtOJ"
      },
      "source": [
        "----\n",
        "### 練習問題\n",
        "\n",
        "次回の授業時にこの問題の答えを解答する小テストを行いますので，事前にやっておいてください．"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "gMFwCD2CebJx"
      },
      "source": [
        "［問題］\n",
        "\n",
        "ほげおくんは，こつこつためたお金で自動車を購入したいと重い，たまに中古車情報を調べている．ほげおくんが気になっている自動車の価格は，調べ始めた日には30万円だったが，それ以降，表のような推移を示している．\n",
        "\n",
        "|$n$|経過週数 $x_n$|価格[万円]|価格差 $y_n$ [万円]|\n",
        "|----:|----:|----:|----:|\n",
        "|1|0| 30 | 0|\n",
        "|2|4|31|1|\n",
        "|3|8|29|-1|\n",
        "|4|20|28|-2|\n",
        "\n",
        "$(x_1, y_1)$ から $(x_4, y_4)$ までの $N=4$ 個のデータについて，経過週数を説明変数$x$，価格**差**を被説明変数$y$として，\n",
        "$$\n",
        "y = ax+b\n",
        "$$\n",
        "という直線を当てはめる回帰分析を実行し，得られた回帰直線によって，ほげおくんが気になっている自動車の価格推移を予測してみよう．\n",
        "\n",
        "以下，簡単のため，回帰係数の推定値の $\\hat{}$ を省略して $a, b$ と表記することにします．"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xiYubW4yhAaF"
      },
      "source": [
        "(1) 上記のデータについて，次の値をそれぞれ求めなさい．\n",
        "\n",
        "$$\n",
        "\\sum_{n=1}^{N}x_n,\\quad\\sum_{n=1}^{N}y_n,\\quad\\sum_{n=1}^{N}x_n^2,\\quad\\sum_{n=1}^{N}x_ny_n\n",
        "$$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "q5nJYX76iGE2"
      },
      "source": [
        "(2) 回帰係数の推定値は，次の正規方程式を満たす．この方程式を解いて，解 $(a,b)$ を求めなさい．\n",
        "\n",
        "$$\n",
        "\\begin{aligned}\n",
        "a\\sum_{n=1}^{N}x_n^2 + b\\sum_{n=1}^{N}x_n &= \\sum_{n=1}^{N}x_ny_n \\\\\n",
        "a\\sum_{n=1}^{N}x_n + b\\sum_{n=1}^{N}1 &= \\sum_{n=1}^{N}y_n\n",
        "\\end{aligned}\n",
        "$$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "90oPdHT-jFRn"
      },
      "source": [
        "(3) (2)の結果を利用して，経過週数から価格差を予測する式を求めなさい．それを用いて，ほげおくんが気になっている自動車の価格の，調べ始めから36週後の価格はいくらと予想されるか答えなさい．"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4IvcDVkpj2j_"
      },
      "source": [
        "［補足］\n",
        "\n",
        "回帰係数は，以下の式の値を計算して求めてもよいのですが，$s_x^2$ や $s_{xy}$ を求めるよりも上記の連立方程式を解くほうがきっと楽ですので，ここではそのような解き方へ誘導しています．\n",
        "\n",
        "$$\n",
        "\\hat{a} = \\frac{s_{xy}}{s_x^2}\\qquad \\hat{b} = \\bar{y} - \\hat{a}\\bar{x}\n",
        "$$"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "pxhkYaaIiAhb",
        "vscode": {
          "languageId": "python"
        }
      },
      "outputs": [],
      "source": []
    }
  ],
  "metadata": {
    "colab": {
      "authorship_tag": "ABX9TyNYOT6qaG6YS5lZjmIBX+ll",
      "collapsed_sections": [],
      "include_colab_link": true,
      "name": "ex09notebookA.ipynb",
      "provenance": []
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}