簡単に機械学習に入る-線形回帰実戦


小文|公衆番号小文のデータの旅
前回は線形回帰の理論部分を紹介しましたが、今回はもちろん期待していた実戦編です!次にstasmodelsパケットの最小二乗法、skleranの最小二乗法、ロット勾配降下法、ランダム勾配降下法、小ロットランダム勾配降下法等から線形回帰を実現する.まず、いくつかの重要な数式を思い出します.
損失関数:最小二乗法最適パラメータを求める:勾配降下法最適パラメータを求める:
次に述べるいくつかの線形回帰を実現する方法は、これらの数式に基づいて得られています.
import pandas as pd
import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score
from sklearn.model_selection import train_test_split
import statsmodels.formula.api as smf

まずデータの長さを見てみましょう.これは、3つのフィーチャーと1つのラベルからなるデータセットです.また,統計的に3つの特徴値の分布が異なり,後続のモデリングに一定の影響を及ぼすことが分かった.第一に、勾配の低下速度に影響を与える.第二に、重みの大きさに影響します.例えば、TVという特徴値は、最大値が296で、平均値が147で、radioの5倍以上で、tvの重みがradioよりはるかに大きくなり、極端にsalesの値はTVに依存していると考えられます.このような状況を回避するために、モデルを構築する前に、データを前処理する必要があります.ここではz-scoreの方法でデータを処理します.
data = pd.read_csv('./Desktop/Advertising.csv',sep = ',')
print(data.describe())
            TV       radio   newspaper       sales
count  200.000000  200.000000  200.000000  200.000000
mean   147.042500   23.264000   30.554000   14.022500
std     85.854236   14.846809   21.778621    5.217457
min      0.700000    0.000000    0.300000    1.600000
25%     74.375000    9.975000   12.750000   10.375000
50%    149.750000   22.900000   25.750000   12.900000
75%    218.825000   36.525000   45.100000   17.400000
max    296.400000   49.600000  114.000000   27.000000

#             ,          
train,test = train_test_split(data,test_size = 0.2,shuffle = True,random_state = 0)
train.iloc[:,:-1] = (train.iloc[:,:-1]-train.iloc[:,:-1].mean())/train.iloc[:,:-1].std()

また,処理後の特徴値モデリングを用いて,各特徴に対応するパラメータ,すなわち上述した重みが各特徴値の重要性を表す.そこで,重み値により有用な特徴をスクリーニングすることができ,Lasso回帰はこの原理に基づいて特徴選択を行う(重みが小さい).では,ある特徴自体が他の特徴よりも重要であれば,その重みを予め大きくしてモデリング,すなわち局所重み付け線形回帰を行うことができる.
局所重み付け線形回帰パラメータの最適解:wは重み行列であり、各特徴に重みを追加する.
ここではまず多くの拡張にすぎず,statsmodelsパケットを介して最小二乗法の線形モデルを構築する.
#statsmodels 、     
stats_model = smf.ols('sales~ TV + radio + newspaper',data = train).fit()
print(stats_model.summary())
OLS Regression Results :                           
==============================================================================
Dep. Variable:                  sales   R-squared:                       0.907
Model:                            OLS   Adj. R-squared:                  0.905
Method:                 Least Squares   F-statistic:                     505.4
Date:                Wed, 19 Jun 2019   Prob (F-statistic):           4.23e-80
Time:                        22:41:19   Log-Likelihood:                -297.29
No. Observations:                 160   AIC:                             602.6
Df Residuals:                     156   BIC:                             614.9
Df Model:                           3                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept     14.2175      0.124    114.463      0.000      13.972      14.463
TV             3.7877      0.125     30.212      0.000       3.540       4.035
radio          2.8956      0.132     21.994      0.000       2.636       3.156
newspaper     -0.0596      0.132     -0.451      0.653      -0.321       0.202
==============================================================================
Omnibus:                       13.557   Durbin-Watson:                   2.038
Prob(Omnibus):                  0.001   Jarque-Bera (JB):               15.174
Skew:                          -0.754   Prob(JB):                     0.000507
Kurtosis:                       2.990   Cond. No.                         1.42
==============================================================================

前回提案したように,F検査によりモデルの良し悪し,T検査判断パラメータが顕著であるか,R方判断引数による因変数の解釈度合いを判断できる.上の図から分かるようにstats_モデルのR方は0.907でフィッティングの程度が良いことを示し,またモデルはF検査に合格したが,パラメータはnewspaper以外はT検査に合格したため,特徴newspaperを除去し,モデルを再構築した.
stats_model1 = sfa.ols('sales~ TV + radio',data = train).fit()
print(stats_model1.summary())
OLS Regression Results :                           
==============================================================================
Dep. Variable:                  sales   R-squared:                       0.907
Model:                            OLS   Adj. R-squared:                  0.905
Method:                 Least Squares   F-statistic:                     761.9
Date:                Wed, 19 Jun 2019   Prob (F-statistic):           1.50e-81
Time:                        22:41:35   Log-Likelihood:                -297.40
No. Observations:                 160   AIC:                             600.8
Df Residuals:                     157   BIC:                             610.0
Df Model:                           2                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept     14.2175      0.124    114.754      0.000      13.973      14.462
TV             3.7820      0.124     30.401      0.000       3.536       4.028
radio          2.8766      0.124     23.123      0.000       2.631       3.122
==============================================================================
Omnibus:                       13.633   Durbin-Watson:                   2.040
Prob(Omnibus):                  0.001   Jarque-Bera (JB):               15.256
Skew:                          -0.756   Prob(JB):                     0.000487
Kurtosis:                       3.000   Cond. No.                         1.05
==============================================================================

stats_model 1はすべて検証に合格したので、式はsales=3.7820 tv+2.8766 radio+14.2175と書くことができ、sklearnのlinear_モデル構築モデル.
#     
x = data.iloc[:,:-2]
y = data.iloc[:,-1:]
x = (x-x.mean())/x.std()
x_train,x_test,y_train,y_test = train_test_split(x,y,test_size = 0.2,
                                      shuffle = True,random_state = 0)

#sklearn      
lr = LinearRegression()
lr.fit(x_train,y_train)
result_lr = lr.predict(x_test)
print('r2_score:{}'.format(r2_score(y_test,result_lr)))  #R  
print('coef:{}'.format(lr.coef_))
print('intercept:{}'.format(lr.intercept_))
r2_score:0.8604541663186569
coef:[[3.82192087 2.89820718]]
intercept:[14.03854759]

sklearnが構築したモデル式はsales=3.8279 tv+2.8922 radio+14.0385である.
次に最小二乗法に基づいてパラメータ最適解を求めることにより,最小二乗法解関数をカスタマイズし,
マトリクスの演算にかかわるため、最初はデータセットをマトリクスのフォーマットに変換します.
#       
def ols_linear_model(x_train,x_test,y_train,y_test):
    x_train.insert(0,'b',[1]*len(x_train))      #       , x0  1    
    x_test.insert(0,'b',[1]*len(x_test))
    x_train = np.matrix(x_train)
    y_train = np.matrix(y_train)
    x_test = np.matrix(x_test)

#          ,         
    if np.linalg.det(x_train.T*x_train) == 0:       
        print('    ,   ')
    else:
#      
        weights = np.linalg.inv(x_train.T*x_train)*x_train.T*y_train 
 #  
        y_predict = x_test*weights                    
        print('r2_score:{}'.format(r2_score(y_test,y_predict)))
        print('coef:{}'.format(weights[1:]))
#  x0 1,           
        print('intercept:{}'.format(weights[0]))      

#  
ols_linear_model(x_train,x_test,y_train,y_test)
r2_score:0.860454166318657
coef:[[3.82192087]
 [2.89820718]]
intercept:[[14.03854759]]

手書き最小二乗法により構築されたモデル式はsales=3.8219 tv+2.882 radio+14.0385である.手書き最小二乗法により,最小二乗法の欠点は,xTxが可逆的でない場合,すなわちxが満ランク行列でない場合に無解であり,この場合最小二乗法により最適解を求めることができないことであることが分かった.そこで最適化アルゴリズムである勾配降下法が得られた.
勾配降下法は、パラメータを更新するたびに使用されるデータセットの数によって、ロット勾配降下法、ランダム勾配降下法、小ロット勾配降下法に分けられる.
  • ロット勾配降下法:パラメータを更新するたびに、すべてのデータセット(all)
  • を使用する
  • ランダム勾配降下法:パラメータを更新するたびに、あるデータセット(one)
  • のみを使用する
  • 小ロット勾配降下法:パラメータを更新するたびに、一部のデータセット(one)を使用する
  • #      :
    def gradient_desc(x_train, y_train,x_test,alpha, max_itor):
        x_train = np.array(x_train)
        x_test = np.array(x_test)
        y_train = np.array(y_train).flatten()
        theta = np.zeros(x_train.shape[1])  
        episilon = 1e-8
        iter_count = 0
        loss = 10
    
    #                       
        while loss > episilon and iter_count < max_itor:   
            loss = 0
            iter_count+=1
    #  (         )
            gradient = x_train.T.dot(x_train.dot(theta) - y_train)/ len(y_train)  
            theta = theta - alpha * gradient 
    #    
            loss = np.sum((y_train - np.dot(x_train, theta))**2) / (2*len(y_train)) 
    
        y_predict = x_test.dot(theta)
        print('r2_score:{}'.format(r2_score(y_test,y_predict)))
        print('coef:{}'.format(theta[1:]))
        print('intercept:{}'.format(theta[0]))
    
    #  
    gradient_desc(x_train, y_train,x_test,alpha=0.001, max_itor=10000)
    r2_score:0.8604634817515153
    coef:[3.82203058 2.8981221 ]
    intercept:14.037935836020237

    構築されたモデル式はsales=3.8220 tv+2.8981 radio+14.0379です.手書きロット勾配降下法により,パラメータを更新するたびに訓練セットのすべてのデータを使用する必要があり,訓練セットのデータ量が大きいと演算時間が長いという欠点が分かった.そこで最適化アルゴリズムであるランダム勾配降下が現れた.
    #      :
    def s_gradient_desc(x_train, y_train,x_test,alpha, max_itor):
        x_train = np.array(x_train)
        x_test = np.array(x_test)
        y_train = np.array(y_train).flatten()
        theta = np.zeros(x_train.shape[1])  
        episilon = 1e-8
        iter_count = 0
        loss = 10
    
    #                       :
        while loss > episilon and iter_count < max_itor:   
            loss = 0
            iter_count+=1
            rand_i = np.random.randint(len(x_train))
    #  (         ):
            gradient = x_train[rand_i].T.dot(x_train[rand_i].dot(theta) - y_train[rand_i]) 
            theta = theta - alpha * gradient
    #    :
            loss = np.sum((y_train - np.dot(x_train, theta))**2) / (2*len(y_train))   
    
        y_predict = x_test.dot(theta)
        print('r2_score:{}'.format(r2_score(y_test,y_predict)))
        print('coef:{}'.format(theta[1:]))
        print('intercept:{}'.format(theta[0]))
        print('iter_count:{}'.format(iter_count))
    
    #  
    s_gradient_desc(x_train, y_train,x_test,alpha=0.001, max_itor=10000)
    r2_score:0.8607601654222723
    coef:[3.83573278 2.90238477]
    intercept:14.036801544903055

    構築されたモデル式はsales=3.8357 tv+2.9023 radio+14.0368です.ランダム勾配降下法を構築することによって,その欠点は,パラメータを更新するたびにあるデータを集中的に訓練するだけで,局所最適解が得られる可能性があることであることが分かった.そこで,ロット勾配降下法とランダム勾配降下法の優位性を総合し,小ロット勾配降下法を得た.
    #       :
    def sb_gradient_desc(x_train, y_train,x_test,alpha,num,max_itor):
        x_train = np.array(x_train)
        x_test = np.array(x_test)
        y_train = np.array(y_train).flatten()
        theta = np.zeros(x_train.shape[1])  
        episilon = 1e-8
        iter_count = 0
        loss = 10
    
    #                       :
        while loss > episilon and iter_count < max_itor:  
            loss = 0
            iter_count+=1
            rand_i = np.random.randint(0,len(x_train),num)
    #  (           ):
            gradient = x_train[rand_i].T.dot(x_train[rand_i].dot(theta) - y_train[rand_i])/num 
            theta = theta - alpha * gradient
    #    :
            loss = np.sum((y_train - np.dot(x_train, theta))**2) / (2*len(y_train))   
       
        y_predict = x_test.dot(theta)
        print('r2_score:{}'.format(r2_score(y_test,y_predict)))
        print('coef:{}'.format(theta[1:]))
        print('intercept:{}'.format(theta[0]))
        print('iter_count:{}'.format(iter_count))
    
    #  
    sb_gradient_desc(x_train, y_train,x_test,alpha=0.001,num=20,max_itor=10000)
    r2_score:0.860623250516056
    coef:[3.82871666 2.89894667]
    intercept:14.042705519319549

    構築されたモデル式はsales=3.8287 tv+2.8989 radio+14.0427です.
    まとめてみます.
  • statsmodels :sales = 3.7820tv + 2.8766radio + 14.2175
  • sklearn:sales=3.8279tv+2.8922radio+14.0385
  • 手書き最小二乗法:sales=3.8219 tv+2.882 radio+14.0385
  • ロット勾配低下:sales=3.8220 tv+2.881 radio+14.0379
  • ランダム勾配降下:sales=3.8357 tv+2.9023 radio+14.0368
  • 小ロット勾配低下:sales=3.8287 tv+2.889 radio+14.0427
  • — end —
    小文のデータの旅
    右上の「+注目」を押して最新shareを取得
    好きならシェアorいいね