Hands-On Machine Learning with Scikit-Learn and TensorFlowを読んでいます．
読んだところをまとめ．

Working with Real Data

不動産会社のデータサイエンティストになったと思って，データ分析してみましょうというストーリー．real-world dataを使いたいということで，今回はCalifornia Housing Prices datasetを使います．著者のGitHubから落とせます．(./datasets/housingにあります．housing.csvというファイル．)

Look at the Big Picture

最初のタスクはカリフォルニア州の調査データからカリフォルニアの住宅価格のモデルを作ることです！このデータには人口，所得の中央値，住宅価格の中央値などがカリフォルニアの地区ごとに収録されてます．その地区はUS Census Breauが提供している最小単位だそう．めんどいので"districts"と呼びます．

Frame the Problem

ボスに最初に聞いとかなきゃいけないのは”ビジネスの目的はなんだい？”ってこと．その回答からこれから作るモデルを会社がどう使って利益を生み出すかということがわかる．ボスの答えは，あなたが作ったモデルの出力が別の機械学習システムに入力されて，最終的にはある地区に投資すべきか否かを判断するのに使うってこと．

次にする質問は，"現状はどうなん？"という質問．この質問に対する回答で参照するべきベンチマークと，問題をどう解決するかの参考情報がわかる．現状は地区ごとの住宅価格が専門家の手作業で作られているようだ．
こいつはコストもかかるし時間もかかる．彼らの算出値もあやしい．これが理由で会社はデータを使って住宅価格を出そうとしているようだ．

以上の情報を元に，frame the problemをしてみよう．
-　教師あり学習？教師なし学習？強化学習？
-　分類問題？回帰問題？それ以外の問題？
-　バッチ学習？オンライン学習？
下に回答が書いてあるけど，それを見る前に自分で考えてみよう．

答え：

教師あり学習(訓練データがすでにある)
多変量回帰問題(目的は予測値を出すというタスク)
バッチ学習(データの規模は小さい)

Select a Performance Measure

いろんな評価指標が書いてある．数式の意味は省略．
RMSE(Root Mean Square Error)
$$\mathrm{RMSE}(\mathbf{X},h)=\sqrt{\frac{1}{m}\sum_{i=1}^m\left(h\left(\mathbf{x}^{(i)}\right)-y^
{(i)}\right)^2}$$
MAE(Mean Absolute Error)
$$\mathrm{MAE}(\mathbf{X},h)=\frac{1}{m}\sum_{i=1}^m\left|h\left(\mathbf{x}^{(i)}\right)-y^
{(i)}\right|$$
高次のノルムの方が大きな値から影響を受け，小さな値からの影響は小さくなる点に注意．なので，RMSEはMAEより外れ値(outlier)の影響を受けやすい．ただし，外れ値があまりなく，分布がベルカーブになっていたらRMSEの方が好まれます．

Check the Assumptions

分析．データの前提をリストアップして検証しとくのはいいことだ．これから作るモデルの出力である地区ごとの住宅は次の機械学習システムに入力されるけど，後工程のシステムは価格を"安い"，"普通"，"高い"というカテゴリーに分けてたらどうでしょう．住宅価格の計算精度を完璧にする意味ないよね？モデルも分類問題として扱った方がよさそう．
今回の例では住宅価格をそのまま使ってるそうなので，回帰問題としてモデルを立てていきましょう．

Get the Data

Create the Workspace

インストールの解説はめんどいのでjupyterが使える前提で話を進めます．
著者のGitHubではなく，自分のGitHubでやってる感じで以下は書いていきます．

Download the Data

import os
import pandas as pd

HOUSING_PATH = "datasets/housing"
def load_housing_data(housing_path=HOUSING_PATH):
    csv_path = os.path.join(housing_path,"housing.csv")
    return pd.read_csv(csv_path)

Take a Quick Look at the Data Structure

　head()を使って取得したデータをちょっと覗いてみよう．
なぜかテキストでは画像貼り付けになっているけど，以下のコードを実行． housing = load_housing_data() housing.head()
各行が1レコード=1地区になっていて，それぞれ10個の属性を持ってます．10個は以下の感じ：longitude, latitude, housing_median_age, total_rooms, total_bedrooms, population, households, median_income, median_house_value, ocean_proximity．
info()を使うとデータの詳細が見れます： housing.info()

ocean_proximityだけは量的情報ではないので，どんなラベルが入っているか見ておきましょう．```
housing["ocean_proximity"].value_counts()

これで内容が見れます．
あと，```describe()```で量的情報のまとめを出力できます．```
housing.describe()

Null値は無視されるのでtotal_bedroomsの数が20433になっている． stdは標準偏差，%の値はN%点を表す(たとえば25%は下から数えて25%点に対応する値が表示される)．
扱っているデータの感覚を掴むためには，もうひとつ別に，ヒストグラムを書いて見るといい．以下のコードを実行する．

import matplotlib.pyplot as plt

housing.hist(bins=50, figsize=(20,15))
plt.show()

このコードで下のようなグラフが出力されます．

ヒストグラムからわかること

median incomeがUSDで表示されていない可能性がある(桁が小さすぎる)．調べてみたところ，スケールか収入がスケール化されていて，最大15最小0.5に制限がかけられている．
housing median ageとmedian house valueも最大値が設定されている．後者は分析の出力となる目標変数なので問題だ，作ったモデルの出力値がこの最大値を声なくなってしまう．クライアントチームと話して，これが問題か議論する必要がある．よくある解決策は以下の2つ：
a. 適切なラベルを集めてくる
b. この地区のラベルを消す
これらの変数はスケールが全く異なっている．
多くのヒストグラムで分布の裾が厚くなっている(medianよりも右側にサンプルが多くなっている)．あとでベル型の分布に変換する方法を考える．
[test setを作る] この段階で訓練データとテストデータを分けておきます．作り方：ランダムにいくつかのインスタンスを抽出して(今回は20%)，訓練データから除けておく．

import numpy as np

def split_train_test(data, test_ratio):
    shuffled_indices = np.random.permutation(len(data))
    test_set_size = int(len(data) * test_ratio)
    test_indices = shuffled_indices[:test_set_size]
    train_indices = shuffled_indices[test_set_size:]
    return data.iloc[train_indices], data.iloc[test_indices]

train_set , test_set = split_train_test(housing, 0.2)
print(len(train_set), " train + ", len(test_set), "test")

これでもいいけど問題あり． →もう一回実行するとテストデータが変わって(ランダムなので)，前回とは別の結果ができてしまう(再現性がない) そういうときは次のようにするといい．

最初に作った訓練データとテストデータをとっておく
np.random.seed(42)などとしてシードを固定すればOK
ただし上の方法は新しいデータセットを加える時に崩壊する．そこで共通の解決策はテストデータに行くか訓練データに行くか決めるインデックスをつくること．こうすることで何回も計算を回してもテストデータはいつも変わらない(=consistent)ようにすることができる．

import hashlib

def test_set_check(identifier, test_ratio, hash):
    return hash(np.int64(identifier)).digest()[-1] < 256 * test_ratio

def split_train_test_by_id(data, test_ratio, id_column, hash=hashlib.md5):
    ids=data[id_column]
    in_test_set = ids.apply(lambda id_: test_set_check(id_, test_ratio, hash))
    return data.loc[~in_test_set], data.loc[in_test_set]

ただし，housingデータセットはID列をもたないので，行番号をIDに使うことにします．

housing_with_id = housing.reset_index()
train_set, test_set = split_train_test_by_id(housing_with_id, 0.2, "index")

行番号をユニークなIDとして使う場合は，新しいデータを既存のデータセットの最後に追加すること，途中の行を削除しないことを確実にしないといけない．これができないのであれば，変更がなさそうな変数を使ってユニークなIDを作ればいい．たとえば，地区のlatitudeと　longitudeは変わらなそうなので，これらをくっつけてIDを作る．

housing_with_id["id"] = housing["longitude"] * 1000 + housing["latitude"]
train_set, test_set = split_train_test_by_id(housing_with_id, 0.2, "id")

Scikit-Learnはデータをいくつかのsubsetsに分解する関数をいくつか用意している．ここではtrain_test_split を紹介．一番簡単な関数でさっきのsplit_train_testと同じことを，追加のオプションをつけて使える．オプションの一つに random_state パラメータというのがあって，このオプションでさっきのシードを指定できる．もう一つの特徴として，行数が同一のデータセットを複数入力することができて，同じインデックスで入力した複数のデータセットを分割することができる．

from sklearn.model_selection import train_test_split
train_set, test_set = train_test_split(housing, test_size=0.2, random_state=42)

ここまでランダムなサンプルを考えてきたけど，この方法だとサンプリングバイアスがかかってしまう．母集団を,strataと呼ばれる同質なsubgroupに分解する層化抽出(stratified sampling)を使えば，テストデータが母集団全体をうまく表すことができるようになる．
median income がmedian housing pricesを予測するのにとても重要な変数である場合，データ全体の収入の各カテゴリをうまく代表してくれているテストデータを作れると嬉しい！このためには各stratumに十分なサンプルが含まれるように層を分けなければならない．言い換えれば，あまり多くの層に分けてしまってはいけないということになる．

housing["income_cat"] = np.ceil(housing["median_income"] / 1.5)
housing["income_cat"].where(housing["income_cat"] < 5, 5.0, inplace=True)

from sklearn.model_selection import StratifiedShuffleSplit

split = StratifiedShuffleSplit(n_splits=1, test_size=0.2, random_state=42)
for train_index, test_index in split.split(housing, housing["income_cat"]):
    strat_train_set = housing.loc[train_index]
    strat_test_set = housing.loc[test_index]
housing["income_cat"].value_counts() / len(housing)

最後に"income_catを消す．あと，結構な時間をテストデータの作成に費やしたけれども，機械学習プロジェクトでは見過ごされがちだがとても重要な部分なんですよ．

for set_ in (strat_train_set,strat_test_set):
    set_.drop("income_cat", axis=1, inplace=True)

Discover and Visulalize the Data to Gain Insights

訓練データを使ってデータの可視化をしてみよう．

Visualizing Geographical Data

housing = strat_train_set.copy()
housing.plot(kind="scatter", x="longitude", y="latitude")
plt.show()

このコードで下のようなグラフが出力されます．

人間の脳は画像情報からパターンを認識する能力が高いんだけど，ここではいろんなパラメータをいじってみよう！

パラメータs：丸の直径が各サンプル点の人口を表すようにする
パラメータc：色で価格を表すようにする
パラメータcmap：cの色に既定カラーセットを使う

housing.plot(kind="scatter", x="longitude", y="latitude", alpha=0.4, s=housing["population"]/100, label="population", figsize=(10,7), c="median_house_value", cmap=plt.get_cmap("jet"), colorbar=True)
plt.show()
plt.legend()

このコードで下のようなグラフが出力されます．

このプロットを使えば，クラスタリングとかに使えそう．

Looking for Correlation

corr_matrix = housing.corr()
corr_matrix["median_house_value"].sort_values(ascending=False)

pandasを使って相関の散布図を作るとさらに便利

from pandas.tools.plotting import scatter_matrix
attributes = ["median_house_value", "median_income", "total_rooms", "housing_median_age"]
scatter_matrix(housing[attributes], figsize=(12,8))
plt.show()

このコードで下のようなグラフが出力されます．

対角成分は相関を出力してもしかたがないので，それぞれの変数のヒストグラムが表示されます． median_house_valueの予測に一番重要な変数はmedian_incomeなのでこいつだけ取り出してみてみよう．

housing.plot(kind="scatter", x="median_income", y="median_house_value", alpha=0.1)
plt.show()

このコードで下のようなグラフが出力されます．

このプロットからわかること
確かに相関は高い(全体的に右上がりの傾向がみれるでしょ？)
価格の最大値の制限がはっきりわかる(median_house_value=500kのところで横線になってる．よくみると350kや450kのあたりにも横線っぽいのがでている．→こういうデータ特性をモデルが出力しないようにするために，これらのデータを取り除いておくというのは一つの手です．

Experimenting with Attribute Combinations

これまでのデータ可視化で取り除きたいデータのくせや裾の厚い分布が変数に含まれていることがわかりました．最後の分析前の下処理として，変数を組み合わせてこうした特徴を消しましょう．たとえば，

total_rooms自体が大事なのではなくて，number_of_rooms_per_houseoldが大事なんだ
total_bedroomsが大事なんじゃなくてbedrooms_per_roomが大事なんだ
pupulation自体が大事なんじゃなくてpopulation_per_householdが大事なんだ
などです．これをやってみよう．

housing["rooms_per_household"]=housing["total_rooms"]/housing["households"]
housing["bedrooms_per_rooms"]=housing["total_bedrooms"]/housing["total_rooms"]
housing["population_per_household"]=housing["population"]/housing["households"]

corr_matrix = housing.corr()
corr_matrix["median_house_value"].sort_values(ascending=False)

Hay, not bad!だそうです．bedrooms_per_roomはtotal_roomsやtotal_bedroomsより高いもmedian_house_valueとの相関になっています．
(明らかに部屋数/ベッド数の比率が低い方が高いでしょ？)
同じような理由でtotal_roomsよりrooms_per_householdの方が情報持っているというのも理解できるでしょう．

こういった探索を完全に行う必要はないです．
大事なことは，早く手をつけてインサイトを集めることが合理的で良いプロトタイプを作ることにつながるってこと．
ただしこれは継続的なプロセスであって，一旦プロトタイプを作ってモデルを作ったとしても，その出力を分析してまたデータの探索ステップに戻る必要があります．

Prepare the Data for Machine Learning Algorithms

データを用意する代わりに，それをする関数を書きましょう．以下のようなメリットがあります

どんなデータセットに対しても簡単に再変換できます．
再利用可能な変換関数のライブラリを作ることができます．
システムで使えます．
いくつかの変換を簡単にトライすることができます．
(説明)変数とラベルに同じ変換をする必要はないので，これらを分けましょう．

housing = strat_train_set.drop("median_house_value", axis=1)
housing_labels = strat_train_set["median_house_value"].copy()

Data Cleaning

total_bedroomsに欠損があるので，まずはこれを直します．オプションは以下の3つ

抜けてる地区(1レコード全部)を消します：dropna()
変数自体を消します：drop()
なにかの値を入れます：fillna()
オプション3でmedianなどを採用した場合，　計算した値(median)を保存しておくのを忘れずに．あとでテストデータを使うときとか，新しいデータがきたときに使います．
オプション3を採用するときは，Scikit-Learnには欠損値を扱う便利なクラス：Imputerがあります．
最初にImputerクラスを作り，リプレイスしたい変数をしていします．
(この操作はデータに関係ない)

from sklearn.preprocessing import Imputer
imputer = Imputer(strategy="median")

変数ocean_proximityはテキスト変数なのでmedianを計算できないため，これを外してからImputerに入れます． Imputerに入れるにはfit()を使います．

housing_num = housing.drop("ocean_proximity", axis=1)
imputer.fit(housing_num)

imputerが計算した値(今回はmedian)はstatistics_メソッドでとれます．さっき見た感じtotal_bedroomsだけに欠損値がありましたが，新しく入ったデータの他の変数に欠損があるかもしれないので，全部の変数のmedianを取っておくのが安全です．

imputer.statistics_
housing_num.median().values

imputerで変換したデータセットを使ってモデルを作ります．結果はNumpy配列になるのでpandasのDataFrameとして使うには，また変換する必要があります．

X = imputer.transform(housing_num)
housing_tr = pd.DataFrame(X, columns=housing_num.columns)

Handling Text and Categorical Attributes

ocean_proximity変数は数値じゃないということで先ほど弾きましたが，テキストを数値に変換してみましょう．
Scikit-LearnにはLabelEncoderというクラスがあってこれができます．

from sklearn.preprocessing import LabelEncoder
encoder = LabelEncoder()
housing_cat = housing["ocean_proximity"]
housing_cat_encoded = encoder.fit_transform(housing_cat)
print(housing_cat_encoded)
print(encoder.classes_)

これは良いものですね．classes_メソッドで数値とテキストの対応も確認できます．ただし一つ問題があって，機械学習のアルゴリズムは2つの値が近いということは2つの値が遠いよりも似ているという判断をします．上の数値を使うと，0と4よりも0と1の方が近いサンプルと認識されてしまいますが，そんなことはないです．この現象に対する解決策はone-hot encoding = 1 of K 符号化法を使うことです． ScikitLearnでは OneHotEncoder を使うことで整数値カテゴリ変数をone-hotベクタに変更することができます． fit_transform()は2次元配列，housing_cat_encodedは1次元配列なので配列を整形する必要があることに注意してください．

from sklearn.preprocessing import OneHotEncoder
encoder = OneHotEncoder()
housing_cat_1hot = encoder.fit_transform(housing_cat_encoded.reshape(-1,1))
housing_cat_1hot

出力がSciPy sparce matrixになっていることに注意してください．この型は数千のカテゴリ変数を扱う時には便利ですし，ほとんど2次元配列として使えます．しかし，諸々の事情からNumPyの密行列として扱いたいときは，toarray()関数を使います．

housing_cat_1hot.toarray()

出力がSciPy sparce matrixになっていることに注意してください．この型は数千のカテゴリ変数を扱う時には便利ですし，ほとんど2次元配列として使えます．しかし，諸々の事情からNumPyの密行列として扱いたいときは，toarray()関数を使います．

housing_cat_1hot.toarray()

ちなみに，以上の2つの変換をLabelBinalizerを使うと一発でできます．注意点はこっちの関数はNumpyの密行列として出力されるので，sparce matrixにするにはsparse_output=Trueとオプションを設定します．

from sklearn.preprocessing import LabelBinarizer
encoder = LabelBinarizer()
housing_cat_1hot = encoder.fit_transform(housing_cat)
housing_cat_1hot

Custom Transformers

Scikit-Learnでいくつか使いでのある変換関数を提供してるけど，自分で変換関数を書く必要があるシチュエーションに出会ったときにはどうすればいいでしょうか．
こんなときもfit(),transform(),fit_transform()関数を自分で作るクラスのメソッドに用意するだけでScikit_Learnの変換関数と同様に(pipelineとかで)使うことができます．最後のfit_transform()は基底クラスとしてTransformerMixinを追加すればＯＫです．
それと，BaseEstimatorを基底クラスとして追加するときは(そしてコンストラクタに *argsと **kargs　に入れない)，get_params()とset_params()の2つの関数が使えます．これはhyperparameterを自動的にチューニングするときに便利な関数です．

例として，前述した変数の合成をする変換クラスを作成してみましょう．

from sklearn.base import BaseEstimator, TransformerMixin
rooms_ix, bedrooms_ix, population_ix, household_ix = 3,4,5,6

class CombinedAttributesAdder(BaseEstimator, TransformerMixin):
    def __init__(self, add_bedrooms_per_room=True): # no *args or **kargs
        self.add_bedrooms_per_room = add_bedrooms_per_room
    def fit(self, X, y=None):
        return self # nothing else to do
    def transform(self, X, y=None):
        rooms_per_household = X[:, rooms_ix] / X[:, household_ix]
        population_per_household = X[:, bedrooms_ix] / X[:, household_ix]
        if self.add_bedrooms_per_room:
            bedrooms_per_room = X[:, bedrooms_ix] / X[:, rooms_ix]
            return np.c_[X, rooms_per_household, population_per_household, bedrooms_per_room]
        else:
            return np.c_[X, rooms_per_household, population_per_household]
attr_adder = CombinedAttributesAdder(add_bedrooms_per_room = False)
housing_extra_attribs = attr_adder.transform(housing.values)

この例では，　add_bedrooms_per_roomというハイパパラメータがあります(規定値はTrue)．
このハイパパラメータを使えば，機械学習アルゴリズムでbedrooms_per_roomという変数を追加するべきか判定できるかも．

一般的に，データ準備の段階で100%確かにと思えないような変数を追加するかどうか選択するハイパパラメータを作れます．
データ準備段階で自動化すればするほど，変数のコンビを自動で試すことができます．
そうすれば，思いもかけなかった変数のコンビがより効率的に見つかるかもしれません．

Feature Scaling

機械学習アルゴリズムは，インプットの数値変数がそれぞれ全然違うスケールになっていると，うまく動かない(いくつかの例外を除いて)．
これを調整するのがfeature scalig．以下の2つがよくある方法です．

min-mac scaling
値のレンジが0から1になるように，変数値のリスケールとシフトを行う．やり方は，最小値(min)を引いて，最大幅(max-min)で割る．Scikit-LearnではMinMaxScalerという変換クラスを使ってこのスケーリングができる．このクラスはfeature_rangeというハイパパラメータを持っていて，レンジの長さを0-1にしたくないときはこのパラメータを使える．
standardization
変数値から平均(mean)を引いて，標準偏差(テキストではvarianceと書いてあるがstandard deviationの間違いと思われる)で割る方法．なので変換後の変数は平均0，分散1になる．min-max scalingと違い，変数値が特定のレンジに収まらないのでneural networkなどの0-1の間の値をインプットにしたいようなアルゴリズムには適さない変換．しかし，外れ値に対して影響を受けにくくなるという特徴がある．Scikit-LearnではStandardScalerという変換クラスを使ってできる．

Transformation Pipelines

いろんな変換を正しい手順で行わなければならないんですが，Scikit-LearnのPipelineクラスというのを使って実装の手助けをしてもらえる．

from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler

num_pipeline = Pipeline([
    ('imputer', Imputer(strategy="median")),
    ('attribsadder', CombinedAttributesAdder()),
    ('std_scaler', StandardScaler()),
])

housing_num_tr = num_pipeline.fit_transform(housing_num)

Pipelineのコンストラクタは手順を流れを定義した name/estimator のペアのリストを引数にとる．ただし，最後以外のestimatorは変換じゃないといけない(=fit_transform()メソッドを持つ関数)．名前は二重のアンダースコアが入ってなければなんでもＯＫ！

(ややこしいが)Pipelineのfitが呼び出されると，Pipelineに定義されている手順通りに，まず最初のfit_transform()が呼び出され，その手順のアウトプットが次の手順のインプットに渡され，次の手順のfit_transform()が呼び出され．．．というのを最後の手順になるまで繰り返す．最後の手順はfit()を呼び出して終わる．

Pipelineは最後のestimatorと同じメソッドを提供する．上の例だと最後のestimatorはStandardScaler()なので，pipelineは(StandardScalerと同じ)transform()メソッドを持っていて，このメソッドを使うと全てのデータに一連ん変換を施す．

さらに，PandasのDataFrameにデータをフィードできると便利だなあと思う．しかし，Scikit-LearnにはPandasのDataFrameを扱う関数クラスなどはないので，それ用の変換クラスを以下で作成する．

from sklearn.base import BaseEstimator, TransformerMixin

class DataFrameSelector(BaseEstimator, TransformerMixin):
    def __init__(self, attribute_names):
        self.attribute_names = attribute_names
    def fit(self, X, y=None):
        return self
    def transform(self, X):
        return X[self.attribute_names].values

この関数を使って以下のようにすると，数値変数のみPandasのDataFrameにすることができる．
PipelineをDataFrameSelectorで数値変数に対してのみ変換を行って，PandasのDataFrameにする．
その後にその他の変数をDataFrameSelectorで抽出してLabelBinarizerを使って変換すればＯＫ！

num_attribs = list(housing_num)
cat_attribs = ["ocean_proximity"]

num_pipeline = Pipeline([
    ('selector', DataFrameSelector(num_attribs)),
    ('imputer', Imputer(strategy="median")),
    ('attribs_adder', CombinedAttributesAdder()),
    ('std_scaler', StandardScaler()),
])

ただし，LabelBinarizerはPipelineに適応していないため，ちょっと工夫して， CustomBinarizerを作ってPipelineに使います．この問題に関するGitHubの質問

class CustomBinarizer(BaseEstimator, TransformerMixin):
    def fit(self, X, y=None,**fit_params):
        return self
    def transform(self, X):
        return LabelBinarizer().fit(X).transform(X)
    
cat_pipeline = Pipeline([
    ('selector', DataFrameSelector(cat_attribs)),
    ('label_binarizer', CustomBinarizer()),
])
num_house = num_pipeline.fit_transform(housing)
cat_house = cat_pipeline.fit_transform(housing)

しかもScikit-LearnのFeatureUnionを使えば，上の2つの処理を1つのパイプラインにまとめることできる．下のコードのように変換クラスのリストをFeatureUnionに投げれば，そのtransform()が呼ばれたときに，各変換クラスのtransform()メソッドを実行して，それぞれのアウトプットが出てから一つに繋げて，結果を返す．

from sklearn.pipeline import FeatureUnion

full_pipeline = FeatureUnion(transformer_list = [
    ("num_pipeline", num_pipeline),
    ("cat_pipeline", cat_pipeline),
])
housing_prepared = full_pipeline.fit_transform(housing)
housing_prepared
housing_prepared.shape

Select and Train a Model

データができたのでとうとうモデル選択と訓練ができるよ！

Training and Evaluating on the Training Set

まずは線形回帰してみましょう．

from sklearn.linear_model import LinearRegression

lin_reg = LinearRegression()
lin_reg.fit(housing_prepared, housing_labels)

この線形回帰モデルを使って予測値と実際の値を比べてみよう！ただし，テキストのコードだと変換がうまくいかないので，housing_preparedから直接データを取ります．

some_data_prepared = housing_prepared[:5]
print("Predictions:", lin_reg.predict(some_data_prepared))
some_labels = housing_labels.iloc[:5]
print("Labels:", list(some_labels))

予測値がそんなに精度よくないっぽいんだけど，モデル全体でどうかわからないので，Scikit-Learnのmean_squared_errorを使って平均平方二乗誤差(RMSE)を計算してみましょう．

from sklearn.metrics import mean_squared_error
housing_predictions = lin_reg.predict(housing_prepared)
lin_mse = mean_squared_error(housing_labels, housing_predictions)
lin_rmse = np.sqrt(lin_mse)
lin_rmse

Okay，なにもないよりはいいけど，great!ってかんじじゃないね．という評価です．なぜなら，各地区のmedian_housing_valueは$120,000〜$265,000というレンジになっているのに，誤差が平均して$68,234もあると精度としては厳しいよね．これはモデルのunderfittingによるものです．こうしたことが起こるということは

変数が提供する情報がいい予測をするのに十分なものになっていない
そもそも今回採用しているモデル(線形回帰モデル)が適していない
という問題があると思われます．ということで解決策は
強力なモデルを採用する
より良い変数を使ってモデルの訓練を行う．
モデルの制約を弱める
というものがあります．今回は1.の方向で解決を探っていきましょう．

線形回帰モデルの代わりにDecisionTreeRegressorを使ってみましょう．
このモデルの詳細は6章で解説します．

from sklearn.tree import DecisionTreeRegressor
tree_reg = DecisionTreeRegressor()
tree_reg.fit(housing_prepared, housing_labels)

housing_predictions = tree_reg.predict(housing_prepared)
tree_rmse = np.sqrt(mean_squared_error(housing_labels, housing_predictions))
print(tree_rmse)

wait what!? No error at all? もちろんそんなことはなくて，これは完全にoverfittingしてしまった結果です．本当にoverfittingしてしまったのかということを次の項で確認します．

Better Evaluation Using Cross-Validation

モデルを評価するひとつの方法はtrain_test_splitを使って訓練データをさらに小さな訓練データとバリデーション用データに分割して訓練し直してみることです．この作業は大変面倒なんですが，Scikit-Learnのcross-validation機能を使うとらくちんにできます．以下のコードではK-fold cross-validationをします．K-fold cross validationとは，訓練データをランダムにK(下のコードでは10)個の部分集合(fold)に分けてモデルをK回訓練します．Kこのうち1つを選んでvalidationに使い，残りのK-1個のデータ集合を使って訓練をします．

from sklearn.model_selection import cross_val_score
scores = cross_val_score(tree_reg, housing_prepared, housing_labels, scoring = "neg_mean_squared_error", cv=10)
tree_rmse_scores = np.sqrt(-scores)

def display_scores(scores):
    print("scores:", scores)
    print("Mean:", scores.mean())
    print("Standard deviation:", scores.std())
    
display_scores(tree_rmse_scores)

これをみるとDecision Treeは前の線形回帰とそれほど良くなっているようには見えないよね．それどころか少し悪化しているように見えます．cross-validationはモデルのパフォーマンスを推定するだけでなく，推定の正確性も評価できるようになるというメリットがあります．
同じスコアを線形回帰モデルも計算してみましょう．

lin_scores = cross_val_score(lin_reg, housing_prepared, housing_labels, scoring = "neg_mean_squared_error", cv=10)
lin_rmse_scores = np.sqrt(-lin_scores)
display_scores(lin_rmse_scores)

ということで線形回帰もDecision TreeもイマイチだったところでRandomForestRegressorを試してみましょう．(これは7章で解説します)

from sklearn.ensemble import RandomForestRegressor
forest_reg = RandomForestRegressor()
forest_scores = cross_val_score(forest_reg,  housing_prepared, housing_labels, scoring = "neg_mean_squared_error", cv=10)
forest_rmse_scores = np.sqrt(-forest_scores)
display_scores(forest_rmse_scores)

Wow, this is much better．
ただし，気をつけなければならないのは，訓練データのスコアはバリデーションデータよりもかなり低くなるため，モデルはまだ訓練データにoverfittingしている可能性がたかいということ．解決策としては

モデルの簡単化
モデルに制約をつける
訓練データをもっと集める
というものがある．また，いろいろ解決策を図る前に，いろんなモデルを試してみることから始めた方が時間の節約になります．

Fine_Tune Your model

モデルのshortlistを手に入れたので，fine-tuneしましょう．

Grid Search

Scikit-LearnのGridSearchCVを使うのが便利です．これは，試したいハイパパラメータを渡せば，ハイパパラメータのすべての組み合わせを試して結果を教えてくれる優れモノなんです．

from sklearn.model_selection import GridSearchCV
param_grid = [
    {'n_estimators': [3, 10, 30], 'max_features': [2, 4, 6, 8]},
    {'bootstrap': [False], 'n_estimators': [3, 10], 'max_features': [2, 3, 4]}
]
forest_reg = RandomForestRegressor()
grid_search = GridSearchCV(forest_reg, param_grid, cv=5,scoring='neg_mean_squared_error')
grid_search.fit(housing_prepared, housing_labels)

param_gridはRandomForestRegressorに入れたいパラメータの組み合わせを入れます．最初の辞書でn_estimator：3個 × max_features：4個 = 12個のハイパパラメータを試します．2個めの辞書で2×3=6のハイパパラメータを試し，合計で18個のハイパパラメータの組み合わせを試します．
結果は以下の感じで取り出せます．

grid_search.best_params_
grid_search.best_estimator_

以下の感じで各試行のスコアもみることができます．

cvres = grid_search.cv_results_
for mean_score, params in zip(cvres["mean_test_score"], cvres["params"]):
    print(np.sqrt(-mean_score), params)