ex04kadaiS.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# AProg2020 ex04 課題S\n",
    "\n",
    "https://www-tlab.math.ryukoku.ac.jp/wiki/?AProg/2020/ex04#kadaiS\n",
    "\n",
    "おまけ課題です．\n",
    "分からないところがあれば随時 takataka に相談してください．"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ねらい\n",
    "\n",
    "この課題のねらいは次のようなものです．\n",
    "- Pythonで初歩的な**データ分析**を経験してみる\n",
    "- ネット上のリファレンス等の情報を自分で収集してプログラムを書く経験をする"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step1  データの準備，読み込み，可視化"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. 以下のリンク先のファイルをこの notebook と同じ場所に保存しましょう: \n",
    "https://www-tlab.math.ryukoku.ac.jp/~takataka/course/AProg/gorigori.csv\n",
    "1. VSCode 等で内容を確認しましょう．\n",
    "\n",
    "このファイルの各行には2つの数が `,`（カンマ） で区切られて並んでいることがわかりますね．このようなファイルの形式を **CSV** (Comma-Separated Values/Variables) といい，拡張子`.csv` のついたファイル名がよく使われます．\n",
    "CSVは，Microsoft Excel 等の表計算ソフトのデータを他のソフトウェアやプログラムで扱うような場合によく使われます．\n",
    "\n",
    "Python のプログラムでこのCSVファイルを読み込んで分析してみましょう．\n",
    "CSVを読み込むには，自分でがんばって全部書く，[標準の csv モジュール](https://docs.python.org/ja/3.8/library/csv.html)を使う，という方法もありますが，ここでは，後の分析で科学技術計算モジュール [NumPy](https://www.numpy.org/) を使うので，NumPy の機能を使うことにしましょう．というわけで..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. NumPy で CSV ファイルを読み込むには， `np.loadtxt` を使うのが簡単です． *Numpy*, *loadtxt*, *CSV* といったキーワードで検索して，使い方を調べましょう．以下のリンク先のリファレンスを読み解くのがかっちょいいです: \n",
    "https://docs.scipy.org/doc/numpy/reference/generated/numpy.loadtxt.html\n",
    "みそは，オプション引数 `delimiter` を使うところ．\n",
    "1. 調べた結果をもとに，入手した CSV ファイルを読み込んで NumPy の配列を得るコードを以下のセルに書きましょう．\n",
    "`loadtxt` の戻り値は，`data` という名前の変数に代入するようにしてください．"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### ここに CSV ファイルを読み込むコードを書く\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "5. 正しく読み込めたか確認するために，以下のセルを順次実行してみましょう．"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#  配列の大きさをタプルとして返す．\n",
    "#  data は 50x2の2次元配列（50行2列の行列）であることがわかる．\n",
    "data.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#  要素番号をこのように指定して配列の要素を参照できる\n",
    "print(data[0, 0], data[0, 1])\n",
    "print(data[1, 0], data[1, 1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(data[1])  # こうすると（0から数えて）1行目だけ取り出して表示"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(data[:3])  # こうすると最初の3行だけ取り出して表示"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(data[:, 1])  # こうすると1列目だけ取り出す"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "6. `data` の0列目を `x` という名前の変数に，1列目を `y` という名前の変数に，それぞれ代入するコードを次のセルに書き，このセルを実行してみましょう．"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# 以下に上記の代入文を書く．\n",
    "\n",
    "\n",
    "\n",
    "plt.plot(x, y, \"o\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "7. matplotlib の使い方を調べて，上記のグラフが次のように描かれるように修正しましょう．\n",
    "    - x 軸の範囲を [0, 40] にする\n",
    "    - y 軸の範囲を [-5, 125] にする\n",
    "    - 点の色を赤にする"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step2 回帰分析／最小二乗法やってみる\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**＜問題設定＞**\n",
    "\n",
    "上記のグラフの元となっているデータは，ほげおくんがある店でアイス「ゴリゴリ君」を売るアルバイトを50日間やった際の，1日毎の「気温」（単位は度，記号 $x$ で表すことにします）と「ゴリゴリ君売上数」（単位は個，記号 $y$ で表すことにします）を表しています．データの0列目およびグラフの横軸が気温，1列目および縦軸が売上数です．\n",
    "\n",
    "50日アルバイトをやったほげおくんは，将来を見込まれてこの店の店長をまかされることになりました．大学でデータ分析や機械学習を学んだほげおくんは，その知識を活かして売上を予測してみようと思い立ちました．気温 $x$ から売上数 $y$ を予測できれば，翌日の予想気温からその日の売上数を予測し，ゴリゴリ君を適切な数仕入れることができるだろう，というわけです．\n",
    "「ゴリゴリ売って大儲けほげ〜」\n",
    "\n",
    "ほげおくんに代わって予測式を求めてみましょう．"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**＜回帰分析／最小二乗法＞**\n",
    "\n",
    "上記のような問題を解くための最も簡単な方法の一つが，**回帰分析**または**最小二乗法**と呼ばれる手法です．\n",
    "「ゴリゴリ君」の例で説明すると，気温 $x$ と売上数 $y$ の関係が $y = ax+b$ という式で表せると仮定して，変数（**パラメータ**といいます） $(a, b)$ として最適な値を求める，というものです．\n",
    "\n",
    "「ゴリゴリ」君データは，気温と売上数のペアが $50$ 個ありました． $N = 50$ とおいて，これらのデータを $(x_n, y_n)$ ($n = 1, 2, \\dots , N$) と表すことにします． $x_n$ が $n$ 番目の気温の値，$y_n$ が $n$ 番目の売上数の値です（配列の要素番号と違って1からはじまってることに注意）．\n",
    "このとき，回帰分析／最小二乗法の問題をより具体的に定式化すると，次のようになります： このデータに直線 $y = ax+b$ をあてはめたときの誤差（**二乗誤差**といいます）を\n",
    "\n",
    "$$ E = \\sum_{n=1}^{N} (y_n - (a x_n + b))^2 $$\n",
    "\n",
    "と表すとき，$E$ を最小にするパラメータ $(a, b)$ を求めなさい．"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "「数値計算法」の授業で学んだ覚えのあるひともいると思いますが，この問題の解はちゃんと存在しそれなりに簡単な式で表せるので，プログラムを書いて計算することができます．しかし，ここでは NumPy の力を借りて最適なパラメータ $(a, b)$ をさくっと求めてみましょう．\n",
    "\n",
    "NumPy で最小二乗法の解を求めるには，関数 `np.linalg.lstsq` を使います: https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.lstsq.html\n",
    "\n",
    "この関数を使うには，少し準備が必要です．\n",
    "気温のデータ $x_1, x_2, \\dots , x_N$ は `x` という名前の配列に格納してあるのでした．"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "これらの値と同数の 1 をならべた $N\\times 2$ 行列を作る必要があるのです．まずは，`x` と同じ大きさで 1 がならんだ配列は，次のようにして作れます: cf. https://docs.scipy.org/doc/numpy/reference/generated/numpy.ones_like.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.ones_like(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`np.vstack` という関数を使うと，複数の配列を垂直方向に(vertically)積み重ねる(stack)ことができます: cf. https://docs.scipy.org/doc/numpy/reference/generated/numpy.vstack.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xx = np.vstack([x, np.ones_like(x)])\n",
    "print(xx)\n",
    "print(xx.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "しかし，この `xx` は上記の通り $2 \\times N$ です．そこで， `xx` の**転置行列**を求めてこれを `XX` という変数に入れることにしましょう．"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ここに xx の転置を XX に代入するコードを書く\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "これで準備が整いました．`XX` と `y` を関数 `np.linalg.lstsq` に渡せば，解が求まります．\n",
    "この関数の使い方を調べて，以下に解を求めるコードを書きましょう．\n",
    "https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.lstsq.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ここに最小二乗法の解を求めるコードを書く．パラメータ (a, b) の値を a, b という変数に代入するようにしよう\n",
    "\n",
    "\n",
    "print(a, b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "得られた $(a, b)$ を使い，$x_n$ ($n = 1, 2, \\dots , N$) を $ax+b$ に代入して得られる値（$ax_n +b$) を求めましょう． `x`, `a`, `b` を使い，結果を `y_est` という名前の配列に代入することにしましょう．"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ここに y_est を計算するコードを書く\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "グラフを描いてみましょう．"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(x, y, \"o\", color = \"red\")\n",
    "plt.plot(x, y_est)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step3 「ゴリゴリ売って大儲けほげ〜」"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "予測式を使って気温から売上数を予想してみよう．"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xr = np.arange(5, 45, 5)  #  区間 [5, 45) の値を 5 刻みで作る\n",
    "print(xr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "以下のセルを修正して正しい売上予測を出せるようにしよう．"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for x in xr:\n",
    "    print('気温 {0} 度のときの売上の予測値は {1:.1f} 個ほげ'.format(x, 999))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "逆に，ある売上数を達成するには気温が何度になればよいか計算してみよう．"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for y in [50, 100, 150, 200]:\n",
    "        print('売上が {0} 個になるには気温が {1:.2f} 度だったらいいほげ'.format(y, 999))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step4  おまけのおまけ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "以下は，読んで実行してみるだけ（$D$ はいろいろ変えてみてね）でokです．\n",
    "\n",
    "**＜最小二乗法による多項式のあてはめ＞**\n",
    "\n",
    "$x$ と $y$ の関係が $D$ 次多項式  $ y = f(x; w_0, w_1, \\dots , w_D) = w_0 + w_1x + w_2x^2 + \\dots +\n",
    "w_Dx^D$ で表せると仮定して，この多項式の$(D+1)$個のパラメータ $w_0, w_1, \\dots , w_D$ として最適な値を最小二乗法で推定してみよう．"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = data[:, 0]\n",
    "y = data[: , 1]\n",
    "\n",
    "# データ数\n",
    "N = len(x)  # x.shape[0]\n",
    "# あてはめる多項式の次数\n",
    "D = 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$N$ 個のデータ $\\{ (x_n, y_n) \\}_{n=1}^{N}$ に対して，次式で表される $(D+1)\\times N$ 行列 $X$ を求めます．\n",
    "\\begin{align}\n",
    "X = \\begin{pmatrix}\n",
    "1 & 1 & \\dots & 1 \\\\\n",
    "x_{1} & x_{2} & \\dots & x_{N} \\\\\n",
    "\\vdots & \\vdots & \\vdots & \\vdots \\\\\n",
    "x_{1}^D & x_{2}^D & \\dots & x_{N}^D \n",
    "\\end{pmatrix}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# X を計算\n",
    "X = np.empty((D+1, N))\n",
    "X[0, :] = 1\n",
    "for d in range(1, D+1):\n",
    "    X[d, :] = x**d\n",
    "print(X.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`X`の転置と `y` を `np.linalg.lstsq` に渡してパラメータを求める．"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rv = np.linalg.lstsq(X.T, y, rcond=None)\n",
    "w = rv[0]  # D+1個のパラメータ\n",
    "print(w)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "関数 $f(x; w_0, w_1, \\dots , w_D)$ をつくる．"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_est = np.poly1d(w[::-1])\n",
    "print(f_est)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "グラフを描く"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xr = np.linspace(5, 35, 100)  # 区間 [5, 35) から等間隔に 100 点をとる\n",
    "\n",
    "plt.plot(x, y, \"o\", color = \"red\")\n",
    "plt.plot(xr, f_est(xr))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "多項式の次数 $D$ をいろいろ変えて実行してみよう．"
   ]
  }
 ],
 "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.8.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {
    "height": "146px",
    "width": "258px"
   },
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "目次",
   "title_sidebar": "Contents",
   "toc_cell": true,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}