機械学習基礎の線形モデル


回帰に使用される線形モデル
せんけいかいき
線形回帰とは何ですか.
回帰問題の一般的なモデルは以下の通りである:$$y=sum w[i]*x[i]+b$$下図に示すように、線形回帰は、与えられた点$(x_i,y_i)$に基づいて直線$$y=ax+b$$をフィッティングし、係数a,bを求める.
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import seaborn as sns
import mglearn
import warnings
warnings.filterwarnings('ignore')

mglearn.plots.plot_linear_regression_wave()
w[0]: 0.393906  b: -0.031804




  は多次元データに広がり,線形回帰モデルの訓練過程はパラメータベクトルwを探す過程であり,フィットしたターゲットが高次元平面になるだけで,線形回帰で最もよく用いられる2つの方法は最小二乗法(OLS)と勾配降下法であり,pythonを用いて線形回帰を実現するにはsklearnとstatsmodelの2つのパケットが利用できる.sklearnは機械学習によく使われるパッケージであり,statsmodelは統計学に偏っている.まず、sklearnのLinearRegressionトレーニングモデルを使用します.
from sklearn.linear_model import LinearRegression 
from sklearn.model_selection import train_test_split 

X,y = mglearn.datasets.make_wave(n_samples=60) #    

#     ,  random_state             
X_train,X_test,y_train,y_test = train_test_split(X,y,random_state=42)
#   sklearn  API,Model().fit(X,y)
lr = LinearRegression().fit(X_train,y_train)

print('  :{}'.format(lr.coef_))
print('  :{}'.format(lr.intercept_))
print('    :{}'.format(lr.score(X_train,y_train)))
print('    :{}'.format(lr.score(X_test,y_test)))
  :[0.39390555]
  :-0.031804343026759746
    :0.6700890315075756
    :0.65933685968637


  sklearnではOLSフィッティングモデルを用い,scoreは可決係数であることから,テストセットの可決係数は0.65程度であり,効果はよくないといえる.これは,元のデータが1次元データであり,データ次元が増加すると線形モデルが非常に強くなるからである.次にstatsmodelを使用してモデルを訓練します.
import statsmodels.api as sm

#        ,     ,          
x = sm.add_constant(X_train)
#    
ols = sm.OLS(y_train,x).fit()
#      
print(ols.summary())
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.670
Model:                            OLS   Adj. R-squared:                  0.662
Method:                 Least Squares   F-statistic:                     87.34
Date:                Sun, 22 Sep 2019   Prob (F-statistic):           6.46e-12
Time:                        21:54:35   Log-Likelihood:                -33.187
No. Observations:                  45   AIC:                             70.37
Df Residuals:                      43   BIC:                             73.99
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const         -0.0318      0.078     -0.407      0.686      -0.190       0.126
x1             0.3939      0.042      9.345      0.000       0.309       0.479
==============================================================================
Omnibus:                        0.703   Durbin-Watson:                   2.369
Prob(Omnibus):                  0.704   Jarque-Bera (JB):                0.740
Skew:                          -0.081   Prob(JB):                        0.691
Kurtosis:                       2.393   Cond. No.                         1.90
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.


 summary出力はeviewsやminitabのような統計スタイルのテーブルであり,可決係数R−squaredは0.67でありsklearn結果と同様であり,モデルのF−statisticおよびパラメータのt値はいずれも結果が顕著であることを示した.
一般最小二乗法(Ordinary Least Square,OLS)
 最小二乗法は、データの実際の値$y_をi$と予測値$hat{y_i}$との偏差が最小、すなわち損失関数が最小であり、OLSは損失関数として平均二乗誤差(Mean Square Error,MSE)を用い、最適化目標は$$minquad MSE=frac 1 nsum_である{i=1}^n(y_i-hat{y_i})^2$$1次元データに対して$$MSE=frac 1 nsum_{i=1}^n(y_i-theta_0-theta_1 x_i)^2$$は最小値を求めるためにバイアスを求める必要がある$frac{partial MSE}{partialtheta_0}=-frac 2 nsum_{i=1}^n(y_i-\theta_0-\theta_1x_i)=0$$,$$\frac {\partial MSE}{\partial\theta_1}=-\frac 2n\sum_{i=1}^n(y_i-\theta_0-\theta_1x_i)x_i=0$$連立可能$theta_1=\frac {\sum(y_i-\bar y)(x_i-\bar x)}{(x_i-\bar x)^2}$$,$$\theta_0=\bar y -\theta_1\bar x$$
ジルコニウムは多重回帰に対して同理であり,以下は行列解法である.損失関数は$J(mathbftheta)=frac{ 1}{2}(mathbf{Xtheta}- mathbf{Y})^T(mathbf{Xtheta}- mathbf{Y})$$導出$$$frac{partial}{partialmathbfmathbftheta} J(mathbftheta)= mathbf{X}^T(mathbf{Xtheta}- mathbf{Y})=$0$最終的に可能である$0$最終的に可能である$0$得$\mathbf{theta}=(mathbf{X^{T}X})^{-1}mathbf{X^{T}Y}$$
  最小二乗原理はMSEを最小化してパラメータ$theta$を求めるために導出することであり,以下では別の方法の勾配降下法を紹介する.
こうばいこうかほう
ジルコニウム勾配降下法は比較的純粋な計算機プログラミング法である.
  図のように、損失関数は係数の関数であることを知っています.一元線形回帰には2つのパラメータがあり、損失関数面を構成しています.まず、上の図面でランダムに初期点を選択し、次の変換$theta_を同時に行う係数のセットをランダムに作成します.0 =\theta_0-\alpha\frac{\partial J(\theta)}{\partial\theta_0}$$$$\theta_1 =\theta_1-alphafrac{partial J(theta)}{partialtheta_1}$$"="が付与番号であり、終了するまでこの手順を繰り返します.
何が起こったのか分析してみましょう.まず,バイアス導関数の係数$alpha$は正数である.バイアス導関数の場合、バイアスがゼロより大きい場合、$J(theta)$は$theta_に従う.i$が大きくなると同時に、新しい$i$は古い$より小さいtheta_i$なので、$J(theta)$が減少します.バイアス導関数がゼロより小さいとき,$J(theta)$は$theta_と共にi$が増加して減少すると同時に、新しい$i$は古い$より大きいtheta_i$,従って,$J(theta)$は減少する,すなわちループごとに損失関数は減少し,最終的には上図に示すように局所的な最小値に達する.
我々の損失関数は凸関数であり,複数の極小値を持つ図形ではなく,その真の図形は以下のように示され,極小値は最小値である.
アルゴリズムの手順:
  • 損失関数
  • を決定する
  • 初期化係数、ステップ長
  • 更新係数
  • 以上の3部を終了
  • まで繰り返す.
    勾配降下法ファミリー
  • バッチ勾配降下法(Batch Gradient Descent)、すなわち前述の方法では、更新のたびにすべてのデータを用いて損失関数および勾配を計算する.
  • ランダム勾配降下法(Stochastic Gradient Descent)は、毎回1つのランダムデータのみを用いて勾配を求める.
  • 小ロット勾配降下法(Mini-batch Gradient Descent)は、一部のデータを用いて勾配を求める.

  • 線形回帰の普及:多項式回帰
      一元線形回帰では,変数yとxが線形関係になっていない場合,線形回帰を直接使用することはできない.テイラーの定理によると:
    a=0で分かるように、yは$x,x^2,x^3...$線形表現は、したがって、$x^n$を追加の変数と見なし、一元線形回帰を多重線形回帰に変換し、モデルの正確性を増加させることができる.
    一般化線形回帰
    すなわち、対数をとることにより、本来の無線形関係の変数を近似線形関係に変換し、線形回帰$$lnmathbf{Y}=mathbf{Xtheta}$$を適用する
    嶺回帰(Ridge)
     嶺回帰はL 2正規化処理回帰モデルを用い、その罰則項はL 2範数であり、罰則項係数は正数であり、sklearnに対応する.Ridgeのパラメータalphaは、alphaを大きくすると係数が0になる傾向があり、トレーニングセットの性能を低下させることは、オーバーフィットを解決する方法であり、同様にsklearn.ridgeはOLSを使用します.  $$J(\mathbf\theta) =\frac{1}{2}(\mathbf{X\theta} -\mathbf{Y})^T(\mathbf{X\theta} -\mathbf{Y}) +\frac{1}{2}\alpha||\theta||_2^2$$
    #    sklearn    
    from sklearn.linear_model import Ridge
    
    X,y = mglearn.datasets.load_extended_boston()
    print('    :{}'.format(X.shape))
    X_train,X_test,y_train,y_test=train_test_split(X,y,random_state = 0)
    
    ridge = Ridge(alpha=1).fit(X_train,y_train)
    LR = LinearRegression().fit(X_train,y_train)
    
    print('      ([  ,  ]):{}'.format([LR.score(X_train,y_train),LR.score(X_test,y_test)]))
    print('     ([  ,  ]):{}'.format([ridge.score(X_train,y_train),ridge.score(X_test,y_test)]))
        :(506, 104)
          ([  ,  ]):[0.9520519609032729, 0.6074721959665752]
         ([  ,  ]):[0.885796658517094, 0.7527683481744755]
    
    

      はbostonデータセットが104個の特徴を有するため、506個のデータしかないことを示す.線形回帰は非常に明らかなオーバーフィットを有し,嶺回帰モデルの訓練精度は線形回帰より低いが,試験精度は線形回帰より高い.
      はまた、データ量を増加させることによってオーバーフィット問題を解決することができ、下図に示すように、データ量が増加すると、線形回帰の試験精度はridgeと類似する.
    mglearn.plots.plot_ridge_n_samples()

    lasso
      lassoはL 1正規化を用い,ペナルティ項はL 1ノルムであるが,ある特徴係数を0にすることができ,モデルをより解釈しやすくすることができ,モデルの重要な特徴を示すこともでき,絶対値を用いることで不導点が存在するためOLS,勾配降下はいずれも利用できない.J(\mathbf\theta) =\frac{1}{2n}(\mathbf{X\theta} -\mathbf{Y})^T(\mathbf{X\theta} -\mathbf{Y}) +\alpha||\theta||_1$$
    #lasso  
    
    from sklearn.linear_model import Lasso
    
    lasso = Lasso(alpha=0.1,max_iter=1000000).fit(X_train,y_train)
    print('    :{}'.format(lasso.score(X_train,y_train)))
    print('    :{}'.format(lasso.score(X_test,y_test)))
    print('       :{}'.format(np.sum(lasso.coef_ !=0)))
        :0.7709955157630054
        :0.6302009976110041
           :8
    
    

    ElasticNet
    ElasticNetはL 1,L 2ノルムを併用して正規化され,$J(mathbftheta)=frac{1}{2 m}(mathbf{Xtheta}-mathbf{Y})^T(mathbf{Xtheta}-mathbf{Y})+alpharho|theta||1 +\frac{\alpha(1-\rho)}{2}||\theta||_2^2$$
    二分類のための線形モデル
    二分類のための線形モデルは、$$y=sum w[i]*x[i]+b>0$$の式で予測できる.
      でよく用いられる二分類モデルには,Logistic回帰(Logistic Regression)とリニアサポートベクトルマシン(Linear Support Vector Machine,リニアSVM)がある.
    Logistic Regression
      Logistic Regressionは,線形回帰により得られる因数yを非線形(Sigmoid関数)変換して[0,1]内にマッピングし,分類サンプル点として0,1の2種類に分類する確率である.
    論理回帰の理解
    Sigmoid関数$$g(z)=frac{1}{1+e^{-z}$$の画像は以下の通りであり,x=0で関数値が0.5,xが無限になる傾向にある場合,関数値はそれぞれ0と1に向かう.$${z=xtheta}$$では線形回帰で得られた関数値を0-1の間にマッピングし,$g(z)$は1に分類される確率と見なすことができ,1に近づくほど1に分類される確率が大きくなり,臨界値0.5で最も誤分類されやすい.
    X = np.linspace(-10,10)
    y = []
    for i in X:
        y.append(1/(1+np.exp(-i)))
    plt.plot(X,y)
    plt.xlabel('z')
    plt.ylabel('g(z)')

    論理回帰の原理
      各サンプルポイントについて$(x_i,y_i)$,$y_i=1,y_i=0$の確率はそれぞれ$$P(y_i=1|x_i,theta)=h_\theta(x_i)$$$$P(y=0|x_i,\theta)=1-h_theta(x_i)$$を$P(y_i|x_i,theta)=h_にマージtheta(x_i)^{y_i}(1-h_theta(x_i)^{1-y_i}$$各サンプル点が独立して同分布すると仮定し、サンプル数はnであり、最大尤度法(MLE)によって尤度関数を構築して$$L(theta)=prod_{i=1}^nP(y_i|x_i,theta)$$尤度関数は既存のサンプルを取得する確率を表しているので最大化すべきであるため,尤度関数の対数の逆数を損失関数$$J(theta)=-lnL(theta)=-sumlimits_{i=1}^{m}(y_iln(h_{\theta}(x_i))+ (1-y_i)ln(1-h_{\theta}(x_i)))$$
    バイアスを求めて$\frac{partial}{partialtheta}J(theta)=X^T(h_{theta}(X)-Y)$$
    勾配降下法$theta=theta-alpha X^T(h_{theta}(X)-Y)$$
    sklearn実装
    #        Logistic Regression
    from sklearn.datasets import load_breast_cancer
    from sklearn.linear_model import LogisticRegression
    
    cancer = load_breast_cancer()
    X_train,X_test,y_train,y_test = train_test_split(cancer.data,cancer.target,stratify = cancer.target,random_state=42)
    
    for C,maker in zip([0.001,1,100],['o','^','v']):
        logistic = LogisticRegression(C = C,penalty='l2').fit(X_train,y_train)
        print('    (C={}):{}'.format(C,logistic.score(X_train,y_train)))
        print('    (C={}):{}'.format(C,logistic.score(X_test,y_test)))
        plt.plot(logistic.coef_.T,maker,label = 'C={}'.format(C))
    plt.xticks(range(cancer.data.shape[1]),cancer.feature_names,rotation = 90)
    plt.xlabel('Coefficient Index')
    plt.ylabel('Coefficient')
    plt.legend()
        (C=0.001):0.9225352112676056
        (C=0.001):0.9370629370629371
        (C=1):0.9553990610328639
        (C=1):0.958041958041958
        (C=100):0.971830985915493
        (C=100):0.965034965034965
    
    
    
    

     Logistic Regressionも正規化を用いることができ,同様に損失関数の後に正規化項を増加させることができる.sklearnではLogisticRegressionはデフォルトでL 2正規化が使用され、パラメータpenaltyは正規化方式を変更できます.上図は異なる正規化パラメータ訓練で得られたモデル係数であり,skleranにおける正規化項Cが小さいほど正規化の程度が強く,パラメータの変換範囲が小さいことがわかる.これは損失関数をとるときに逆数を取ったためであり,−Cはlasso中の$alpha$,Cが小さいほど−Cが大きくなる.$$J(\mathbf\theta) = -\sum\limits_{i=1}^{m}(y_iln(h_{\theta}(x_i))+ (1-y_i)ln(1-h_{\theta}(x_i))) -\frac{1}{2}C||\theta||_2^2$$
    リニアSVC
    SVM参照
    マルチカテゴリ用の線形モデル
     多くの線形モデルは多分類問題には適用されず、データがA、B、Cの3種類に分類されるように、1対の残りの方法を使用することができ、3つの分類器がそれぞれ3つの分類に対応するように訓練する必要がある.例えば、Aクラスの分類器はデータをAクラスとAクラスではないに分け、同時に複数の分類(AクラスとBクラスに属するように)に属し、いずれのクラスにも属しないデータに対して、得点の高いカテゴリに分けられます.
    from sklearn.datasets import make_blobs
    
    X,y=make_blobs(random_state=42)
    mglearn.discrete_scatter(X[:,0],X[:,1],y)
    plt.xlabel('Feature 0')
    plt.ylabel('Featrue 1')
    plt.legend(['Class0','Class1','Class2'])

      上記データに対して線形SVMを用いて分類する:
    from sklearn.svm import LinearSVC
    LSVC=LinearSVC().fit(X,y)
    LSVC.coef_,LSVC.intercept_,LSVC.score(X,y)
    (array([[-0.17492558,  0.23141285],
            [ 0.4762191 , -0.06937294],
            [-0.18914557, -0.20399693]]),
     array([-1.0774515 ,  0.13140521, -0.08604887]),
     1.0)
    
    
    

    LinearSVCは3つの直線を出力し、各直線は1つのカテゴリを他のカテゴリと分離し、以下は可視化し、3つの直線は7つの領域に分け、交差領域は平均的に分配される.
    mglearn.plots.plot_2d_classification(LSVC,X,fill=True,alpha=.7)
    mglearn.discrete_scatter(X[:,0],X[:,1],y)
    line = np.linspace(-10,10)
    for coef,intercept,color in zip(LSVC.coef_,LSVC.intercept_,['b','r','g']):
        plt.plot(line,-(line*coef[0]+intercept)/coef[1])
    plt.xlabel('Feature 0')
    plt.ylabel('Featrue 1')
    plt.legend(['Class0','Class1','Class2','Line class 0','Line class 1','Line class 2'],loc=(1.01,0))