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