<jupyter_start><jupyter_text>- 最小二乗法(ordinary least squares: OLS)による回帰をNumpy, sklearnで．  - 特徴量の方がサンプルデータより多いとOLSできないので，lasso, ridge, elastic netsという3つの手法を試す．- レコメンドという，応用範囲が広い重要な分野の基礎を学ぶ．  まとめ- 最小二乗法の復習    - sklearnのボストン物件価格データを利用    - 特徴量一つだけよりバイアスを加えた方が，さらにたくさん特徴量があった方が一般的に良い精度が出る．    - 回帰分析でも交差検定をして正しい汎化能力を確かめよう．    - 特徴量の数P > サンプル数N のとき，訓練誤差をほぼゼロにすることができるが，汎化能力が著しく落ちる．        - 10-K reportsのデータはP>Nとなっており，これでP>N問題を確かめた．    - 過学習対策        - Lasso: L1罰則項(回帰係数の絶対値の総和)をペナルティとする．            - スパースなモデルに．        - Ridge: L2罰則項(回帰係数の2乗の総和)をペナルティとする．        - Elastic Net: L1, L2の比率を設定し，それらの和を罰則項とする．             - P>N問題の解消，相関の高い特徴量からの複数考慮といった効果がある．- パラメータ探索には2段階の交差検定が必要    - 1段目： テストデータと訓練データに分ける    - 2段目： 1段目の訓練データをさらに分割し，それぞれ違ったパラメータを与えて訓練した結果，誤差が最少となったパラメータを採用- 映画の推薦問題を罰則付き回帰分析で解いてみる    - Netflix Challengeのアカデミック用データを利用    - ターゲットユーザ以外のユーザが行なった評価を特徴量としてLassoで解いてみた    - 全体の平均をスコアの予測値とするよりはLassoを使った方が8％ほどの精度改善が見られた．  回帰を用いて物件価格を予想するボストンの物件の価格を予想する問題を考える．  sklearnのデータセットを使って教師あり学習を行う．<jupyter_code>from sklearn.datasets import load_boston
boston = load_boston()
print(boston.DESCR)
print(boston.feature_names)<jupyter_output>.. _boston_dataset:

Boston house prices dataset
---------------------------

**Data Set Characteristics:**  

    :Number of Instances: 506 

    :Number of Attributes: 13 numeric/categorical predictive. Median Value (attribute 14) is usually the target.

    :Attribute Information (in order):
        - CRIM     per capita crime rate by town
        - ZN       proportion of residential land zoned for lots over 25,000 sq.ft.
        - INDUS    proportion of non-retail business acres per town
        - CHAS     Charles River dummy variable (= 1 if tract bounds river; 0 otherwise)
        - NOX      nitric oxides concentration (parts per 10 million)
        - RM       average number of rooms per dwelling
        - AGE      proportion of owner-occupied units built prior to 1940
        - DIS      weighted distances to five Boston employment centres
        - RAD      index of accessibility to radial highways
        - TAX      full-value property-tax rate per $10,000
        - PTRATIO  pu[...]<jupyter_text>まずは部屋の数の平均値と家の値段をプロットしていく．<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

x = boston.data[:, 5]

plt.scatter(x, boston.target)

x = np.array([ [v] for v in x])
y = boston.target
slope, _, _, _ = np.linalg.lstsq(x, y, rcond=None) #[傾き, バイアス(あれば)], 誤差, ランク, 特異値. rcondはワーニング対策

ran = np.arange(0, 10, 0.1)
plt.scatter(ran, slope * ran)<jupyter_output><empty_output><jupyter_text>あまりうまくフィッティングできていない．  特徴量の要素が1つだけなので，家の値段 × 傾き = 家の値段 という式を最適化していることになるが，これは雑すぎる．  次に，バイアス項を追加してみる．<jupyter_code>x = boston.data[:, 5]

plt.scatter(x, boston.target)

x = np.array([ [v, 1] for v in x]) # バイアス項を追加するため，1を追加 y = ax + b = a x + b 1, この xと1に当たるところをlstsqに渡す．
y = boston.target
(slope, bias), total_error, _, _ = np.linalg.lstsq(x, y, rcond=None)

ran = np.arange(0, 10, 0.05)
plt.scatter(ran, slope * ran + bias)

plt.xlim(3, 10)
plt.ylim(0, 55)

rmse = np.sqrt(total_error[0]/len(x))
print(rmse)<jupyter_output>6.6030713892225625<jupyter_text>外れ値(outlier)にはうまく対応できていないが，結果は良くなっている．  このモデルがどれだけうまくデータに適合できているか，定量的な評価を行うには，平均二乗平方根誤差(root mean squared error: RMSE)を求める．  バイアス項 あり：なし = 7.6:6.6 と改善されている．  RMSEは標準偏差のようなものなので，誤差が正規分布に従うと仮定すると，ほとんどのデータはおおよそ平均値から標準誤差2つ分の範囲に収まる．  -2RMSE ~ 2RMSEの範囲は信頼区間(confidence interval)に相当する．  そのため，この結果からわかることは，我々が予測する家の値段の誤差は多くても$\pm 6.6 \times 2 \times1000 = \pm 13000$ドルの範囲に収まるであろうと言える．  多次元回帰データを全て使って予測してみる．<jupyter_code>x = np.array([ np.concatenate((v, [1])) for v in boston.data])
y = boston.target
slope, total_error, _, _ = np.linalg.lstsq(x, y, rcond=None)
rmse = np.sqrt(total_error[0]/len(x))
print(rmse)<jupyter_output>4.679191295697282<jupyter_text>データをたくさん使ったためか，誤差が改善された．  ここでは14次元の回帰を行なっているため，結果のプロットはできない．  回帰における交差検定回帰における訓練誤差とテスト誤差の差はあまり大きくないことが予想されるが，それでも交差検定を行う．  ここからは，今後より進んだ手法への置き換えをスムーズにするために，sklearnのLinearRegressionを利用する．<jupyter_code>from sklearn.linear_model import LinearRegression

def rmse(model, data, target):
    p = model.predict(data).ravel() # 予測にはmapを使わず，直接データ全体を入力するようになった
    e = p - target # インスタンスごとの誤差を得る
    total_error = np.sum(e * e)
    rmse = np.sqrt(total_error / len(p))
    return rmse
lr = LinearRegression(fit_intercept=True) # fit_intercept: バイアス項を追加
lr.fit(x, y)
print(f"RMSE of training: {rmse(lr, x, y)}")<jupyter_output>RMSE of training: 6.603071389222561<jupyter_text>numpyのRMSEとだいたい同じになった.  KFoldによって10分割の交差検定を行う．<jupyter_code>from sklearn.model_selection import KFold

def kfold_regression(model, x, y, n_splits=10):
    kf = KFold(n_splits=n_splits)
    rmses = []
    for train, test in kf.split(x):
        model.fit(x[train], y[train])
        rmses.append(rmse(model, x[test], y[test]))
    return np.asarray(rmses).mean()

print(f"RMSE on 10-fold CV: {kfold_regression(lr, x, y)}")<jupyter_output>RMSE on 10-fold CV: 5.941570588762131<jupyter_text>訓練誤差よりもテスト誤差は大きくなったが，汎化能力についてより正しい評価を行なっていると言える．  最小二乗法はモデルが単純で高速なため，最初に試してみるモデルとして良い．  罰則付き回帰通常の回帰では過学習を起こすので，パラメータが過剰に適合することに対して罰則を追加する罰則付き回帰(penalized regression)を行なってみる．  罰則付き回帰はバイアス-バリアンスのトレードオフの一つとみることができる．  罰則を用いることは誤差(バイアス)を大きくすることであり，それはバリアンスを減らすことになるので訓練データに適合しすぎることを防ぐ．  よって全体の汎化能力はより高くなると期待できる．    L1, L2罰則項まず，最小二乗法の最適化の式は以下$$ \vec{b^*} = \arg \min_{\vec{b}} (y - X \vec{b})^2 $$L1罰則項(Lasso回帰): 係数の絶対値の和を追加， 誤差を小さくするだけでなく，最小限の係数で抑えようとする．->スパースなモデルに!(不要特徴量削減)$$ \vec{b^*} = \arg \min_{\vec{b}} (y - X \vec{b})^2 + \lambda \sum_i | b_i | $$L2罰則項(Ridge回帰): 係数の二乗の和を追加$$ \vec{b^*} = \arg \min_{\vec{b}} (y - X \vec{b})^2 + \lambda \sum_i b_i^2 $$ここで$\lambda$罰則項の強さを調整するパラメータ．  Lasso回帰とRidge回帰を合わせるとElastic netと呼ばれる手法になる．  $$ \vec{b^*} = \arg \min_{\vec{b}} (y - X \vec{b})^2 + \lambda \sum_i \Bigl\{ \alpha |b_i| + (1-\alpha)b_i^2 \Bigr\} $$詳しくはこちら http://tekenuko.hatenablog.com/entry/2017/11/18/214317  scikit-learnのLassoとElastic netを使用するsklearnで先ほどの最小二乗法をElastic netによる回帰に簡単に差し替えられる．  L1罰則項を追加するにはLasso, L2罰則項を追加するにはRidgeクラスを使えばいい．<jupyter_code>from sklearn.linear_model import ElasticNet
en = ElasticNet(fit_intercept=True, alpha=0.5)
print(f"RMSE on 10-fold CV: {kfold_regression(en, x, y)}")<jupyter_output>RMSE on 10-fold CV: 6.596253559088761<jupyter_text>RMSEは5.18から5.11に減少．    PがNより大きい場合特徴量の数Pがサンプル数Nより大きい場合の問題について考える．  例えば，テキストデータの単語1つ1つを特徴量とすると，特徴量の方が多くなりがちである．  サンプルより特徴量の方が多い時，訓練データに対して完全に当てはめることが可能になり，訓練誤差を0にできる．  なぜならば，連立方程式の数が変数の数より地位し場合，常に連立方程式の解を求められることと同じだからである．  しかしそれは，汎化能力が低くなってしまう原因となる．  テキストに基づいたサンプルデータカーネギーメロン大学のNoah Smith教授らは米証券取引委員会(America SEC)に提出される決算報告書，10-K reportsのデータマイニングを研究している．  10-K reportsを元に株価の推移を予測したい．  16087のサンプルに対して，単語の種類を特徴量とすると150360個にも及ぶ．  データは次のシェルスクリプトでダウンロードできる．(mac, wgetはbrew install wget) ``` bash!/usr/bin/env bashwget http://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/regression/E2006.train.bz2bunzip2 E2006.train.bz2``` データセットはSVMLightというSVM関連のツール用のフォーマットになっており，sklearnでロード可能<jupyter_code>from sklearn.datasets import load_svmlight_file
data, target = load_svmlight_file('./data/E2006.train') # 10秒くらいかかる<jupyter_output><empty_output><jupyter_text>targetは単純な1次元のベクトルだが，dataはスパースな行列．  まずはtargetの各要素の値を見てみる．<jupyter_code>print(f'Target shape: {target.shape}')
print(f'Min target value: {target.min()}')
print(f'Max target value: {target.max()}')
print(f'Mean target value: {target.mean()}')
print(f'StdDev target: {target.std()}')<jupyter_output>Target shape: (16087,)
Min target value: -7.89957807346873
Max target value: -0.519409526940154
Mean target value: -3.5140531366944456
StdDev target: 0.6322783539114604<jupyter_text>一応dataの形も見ておく<jupyter_code>print(f'Data shape: {data.shape}')
# print(f'Data example: {data[0]}')<jupyter_output>Data shape: (16087, 150360)<jupyter_text>とりあえずこのまま最小二乗法にぶち込んでみる<jupyter_code>lr.fit(data, target)
rmse(lr, data, target)<jupyter_output><empty_output><jupyter_text>誤差はかなり小さいが，これはtargetの標準偏差0.6に比べてとても小さい値になっている．  特徴量の数がサンプル数よりはるかに大きいので，過剰適合を起こしていると考えられる．  試しに，交差検定を行なってみる．<jupyter_code>kfold_regression(lr, data, target)<jupyter_output><empty_output><jupyter_text>誤差は0.78程度になる．  しかしそもそも標準偏差が0.6なので，常に-3.5と予測していてもRMSEは0.6程度にしかならないことを考えれば，より悪い結果と言える．  この過学習への対策として正則化が行われる．  elastic netを用いて交差検定を行ってみる．<jupyter_code>en = ElasticNet(alpha=0.5)
en.fit(data, target) 
kfold_regression(en, data, target)<jupyter_output><empty_output><jupyter_text>常に平均値を予測結果とするよりは，良い結果になった．あんまり良くなってないけど．    ハイパーパラメータの賢い設定方法罰則パラメータは大きい値に設定した場合，極端な例としては全ての係数が0になるのが最適解となってしまうため，未学習の可能性がある．  小さすぎると今度はただの最小二乗法になってしまうので，過学習のおそれがある．  最適なパラメータを求めるには，パラメータについて取りうる候補を用意し，交差検定で最適な値を一つ選び出す．  したがって，訓練データとテストデータに分ける交差検定と，最適なパラメータを選び出す交差検定の2段階の交差検定を用いなければならない．  このとき，パラメータ調整のためのデータ(2段目)は1段目の訓練データのみを使うことに注意する．  2段目で分割されたデータはそれぞれ異なるパラメータで訓練され，もっとも結果が良かったパラメータを使ってテストを行った結果を採用する．  訓練やパラメータ調整で利用されたデータは汚されたデータと呼ばれ，汎化能力を計測するために使うことはできないことに留意する．  こういった2段階の交差検定は，LassoCV, RidgeCV, ElasticNetCVで行うことができる．<jupyter_code># CV設定しろっていうワーニングをとりあえず止める
import warnings
warnings.filterwarnings('ignore')
from sklearn.linear_model import ElasticNetCV

met = ElasticNetCV(alphas=[.125, .25, .5, 1., 2., 4.]) # alphaは探索対象なので探す必要がないが，時間かかりすぎなので制限
kf = KFold(n_splits=10)
err = 0
for n, (train, test) in enumerate(kf.split(data)):
    met.fit(data[train], target[train])
    p = met.predict(data[test]).ravel()
    e = p - target[test]
    e = np.dot(e, e)
    err += e
    print(f"Cross Validation: {n}, err: {e}")
print(np.sqrt(err / len(target)))<jupyter_output>Cross Validation: 0, err: 279.76675152278665
Cross Validation: 1, err: 244.654851577256
Cross Validation: 2, err: 249.28855645149685
Cross Validation: 3, err: 258.30418227923565
Cross Validation: 4, err: 296.1305187299544
Cross Validation: 5, err: 249.7932932467228
Cross Validation: 6, err: 247.76458906671712
Cross Validation: 7, err: 271.0096111522141
Cross Validation: 8, err: 280.27278665256557
Cross Validation: 9, err: 299.01610344895033
0.4078548439101767<jupyter_text>レーティング予測とレコメンド映画のお気に入りを予測する問題が「Netflix Challenge」として有名．  コンテストの勝者は，先進の機械学習を組み合わせた手法＋データの前処理への多くの労力によって優勝した．  前処理の段階で，あるユーザが全体的に高い評価をしていたり，常にネガティブな評価をする傾向がある時を考慮しなければならない．  他にも，どれくらい古い映画か？何人のユーザから評価をつけられたか？といった特徴量を利用していた．  良い結果を出すためには，優れたアルゴリズムだけではなく，「手を汚す」必要がある．  コンテストで実際に使われたデータセットを手に入れることはできない．  しかし，同様なデータがアカデミック用として，GroupLensというミネソタ大学の研究所によって公開されている．  クラス分類問題によるアプローチを取ることもできるが，この問題を回帰問題として考え，本章で学んだ手法を適用してみる．  クラス分類問題だと，次の2点を考慮するのが難しいため，回帰問題として扱う方が理にかなっている．  - 予測の誤差は大きいほど重大に評価しなければならない．例えば，4点と5点の誤差は大きくないが，1点と5点の誤差は甚大である．  - 中間的な値が意味を持つ．ユーザの入力は整数だが，予測値は4.7のような少数になっても意味がある．我々には2つの選択肢がある．  - 映画を対象とするモデル- ユーザを対象とするモデル(今回はこっち)      あるユーザを対象として，そのユーザの映画の評価方法を学習する．予測器はある映画を入力すると，予測スコアを出力する．      映画の特徴量は，対象ユーザ以外の他ユーザの映画に対してのスコアとする．      対象ユーザと評価のつけ方が似ているユーザに対して大きな回帰係数がつくことになる． まずはデータを取得する．シェルスクリプトは以下``` bash!/usr/bin/env bashwget http://www.grouplens.org/system/files/ml-100k.zip --no-check-certificateunzip ml-100k.zip```次にデータをロードし，(ch07/usermodel.pyを参照)<jupyter_code>from scipy import sparse # 2次元スパース行列のためのパッケージ

data = np.array([ [int(tok) for tok in line.split('\t')[:3]] for line in open('./data/ml-100k/u.data')])
ij = data[:, :2]
ij -= 1  # original data is in 1-based system
values = data[:, 2]
reviews = sparse.csc_matrix((values, ij.T)).astype(float) # Compressed Sparse Column matrix
print(reviews.shape)
# print(reviews[0])<jupyter_output>(943, 1682)<jupyter_text>学習するデータは疎行列で，評価付けされた映画に対してのみ1~5の値がついている．それ以外は0.  今回は前と同じではつまらないのでLassoCVクラスを利用してみる．<jupyter_code>from sklearn.linear_model import LassoCV
reg = LassoCV(alphas=[.125, .25, .5, 1., 2., 4.]) # 交差検定で用いるalphaの値を制限<jupyter_output><empty_output><jupyter_text>インデックスがi番目のユーザuに対して，モデルを学習する関数を書いていく．  uおための映画の評価方法をモデル化するには，uが評価付けした映画に対してだけ興味を持てば良い．  uによって評価付けされた映画のインデックスを抜き出す．<jupyter_code>def target_user_data(user_id):
    i = user_id
    u = reviews[i]
    u = u.toarray().ravel() # スパース行列を通常の配列に戻してflatten
    ps, = np.where(u > 0) # 評価付けされたインデックスのみを取得
    us = np.delete(np.arange(reviews.shape[0]), i) # ユーザi以外のユーザ配列を作成
    x = reviews[us][:, ps].T # ユーザi以外のユーザがつけた映画のスコア
    y = u[ps] # 対象ユーザuが評価付けを行なった映画のスコアを抜き出す
    return x, y
# 優れた映画は平均値が高くなるため，映画ごとに正規化(平均値を引く) を行う
def movie_norm(xc): # xc: 全ての映画の特徴量(ユーザi以外の評価)
    xc = xc.copy().toarray() # スパースでない配列に変換することでsklearnが高速に処理できる．(メモリ領域の許す限りは)
    x1 = np.array([xi[xi > 0].mean() for xi in xc]) # 全ての映画について，評価の平均値をとる． 評価がついているもののみを対象とする．
    x1 = np.nan_to_num(x1) # 評価が1つもついておらず，NANになる場合があるので，これを0に置き換える．

    # 全ての映画に対し，それぞれの映画の評価の平均値を引く．
    for i in range(xc.shape[0]):
        xc[i] -= (xc[i] > 0) * x1[i] # 平均値を引くが，評価付けしていない映画のスコアは0のまま
        
    return xc, x1 # 正規化した配列, 平均値
def learn_for(user_id, model=reg):
    x, y = target_user_data(user_id)
    kf = KFold(n_splits=10)
    err = 0
    for n, (train, test) in enumerate(kf.split(x)):
        
        xc, x1 = movie_norm(x[train])
        model.fit(xc, y[train] - x1) # 学習にはラベルの評価も正規化を行う．
        
        xc, x1 = movie_norm(x[test])
        p = model.predict(xc).ravel()
        
        e = np.sum( ((p + x1) - y[test]) ** 2 )
        err += e
        # print(f"Cross Validation: {n}, err: {e}")
        
    return np.sqrt(err / len(target))
rmses = []
for i in range(reviews.shape[0]):
    rmse = learn_for(i)
    # print(f"RMSE for user {i}: {rmse}")
    rmses.append(rmse)

print(f"Average RMSE: {np.asarray(rmses).mean()}")<jupyter_output>Average RMSE: 0.06882747921849655<jupyter_text>ユーザのある映画への評価を予測する時，もっとも単純なのは平均値を予測すること．  そのような単純な手法よりは80%良い結果になる．  ここではユーザごとに完全に別のモデルを学習しており，その結果は目覚ましいとは言えない．  次回は回帰の手法を上回る，全てのユーザと全ての映画の情報を統合した賢い方法を紹介する．<jupyter_code># 単純に平均を予測値としてみる
class AverageModel():
    def fit(self, x, y):
        pass
    
    def predict(self, x):
        return np.zeros(x.shape[0])
rmses = []
for i in range(reviews.shape[0]):
    rmse = learn_for(i, model=AverageModel())
    rmses.append(rmse)

print(f"Average RMSE: {np.asarray(rmses).mean()}")
0.06882747 / 0.07407938<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12082636211537062, 'fraction_numerical': 0.048155961300647414, 'mean_word_length': 5.680272108843537, 'pattern_counts': {'":': 0, '<': 62, '<?xml version=': 0, '>': 70, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '18', 'dir': 'jupyter-structured-clean-dedup', 'id': '80514', 'license': 'MIT', 'n_tokens_mistral': 9186, 'n_tokens_neox': 8635, 'n_words': 870, 'path': '7章 レコメンド.ipynb', 'repo_name': 'fuurin/BuildingMachineLearningSystems', 'size': '13747'}