<jupyter_start><jupyter_text>実験Ⅳ－１　正則化付き重回帰モデル1. サンプルデータの読み込み1. 線形重回帰の実行1. 正則化付き基底展開重回帰の実行<jupyter_code>!pip install japanize-matplotlib
%matplotlib inline

import warnings
warnings.simplefilter('ignore') #警告を無視（コメントアウト推奨）

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import japanize_matplotlib #日本語化matplotlib
sns.set(font="IPAexGothic") #日本語フォント設定<jupyter_output>Requirement already satisfied: japanize-matplotlib in /usr/local/lib/python3.6/dist-packages (1.1.2)
Requirement already satisfied: matplotlib in /usr/local/lib/python3.6/dist-packages (from japanize-matplotlib) (3.2.2)
Requirement already satisfied: numpy>=1.11 in /usr/local/lib/python3.6/dist-packages (from matplotlib->japanize-matplotlib) (1.18.5)
Requirement already satisfied: kiwisolver>=1.0.1 in /usr/local/lib/python3.6/dist-packages (from matplotlib->japanize-matplotlib) (1.2.0)
Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.6/dist-packages (from matplotlib->japanize-matplotlib) (0.10.0)
Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.1 in /usr/local/lib/python3.6/dist-packages (from matplotlib->japanize-matplotlib) (2.4.7)
Requirement already satisfied: python-dateutil>=2.1 in /usr/local/lib/python3.6/dist-packages (from matplotlib->japanize-matplotlib) (2.8.1)
Requirement already satisfied: six in /usr/local/lib/python3.6/dist-packa[...]<jupyter_text>Google Colaboratory でファイルを読み込む以下の方法があります。1. ファイルをアップロードする1. URL を指定する<jupyter_code># ファイルをアップロードする
#from google.colab import files
#uploaded_file = files.upload()
#file_name = next(iter(uploaded_file))
#print(file_name)

# URLを指定する
filename = 'https://raw.githubusercontent.com/msiplab/EicEngLabIV/master/data/sample01_02.csv'

#観測データのファイル
#filename = './data/practice01_02.csv'
#filename = './data/sample01_02.csv'
#観測データ（空白区切り）の読み込み
#dataset = pd.read_csv(filename, delim_whitespace=True)
dataset = pd.read_csv(filename)
display(dataset)<jupyter_output><empty_output><jupyter_text>データの整形以下では、二次元テーブルデータの読み込みを想定しています。1. 説明変数の整形（メッシュグリッド化）1. 目的変数の抽出<jupyter_code># 説明変数x1（ヘッダの読み込み）
x1 = dataset.columns[1:].astype(float).values
display(x1)
# 説明変数x2（１列目の読み込み）
x2 = dataset.iloc[:,0].astype(float).values
display(x2)

# 説明変数x1,x2の設定（メッシュグリッド）
x1,x2 = np.meshgrid(x1,x2) 
display(x1)
display(x2)

# 目的変数yの設定
y = dataset.iloc[0:,1:].astype(float).values
display(y)

# 総観測点数の抽出
nSamples = y.size<jupyter_output><empty_output><jupyter_text>観測データの散布図- データの値の範囲を抽出- (x1,y), (x2,y)の散布図を描画<jupyter_code>from mpl_toolkits import mplot3d

# 散布図の範囲情報の抽出
minx1 = np.min(x1)
maxx1 = np.max(x1)
rx1 = maxx1-minx1
cx1 = 0.5*rx1
minx2 = np.min(x2)
maxx2 = np.max(x2)
rx2 = maxx2-minx2
cx2 = 0.5*rx2

# 散布図の表示
fig, axs = plt.subplots(1, 2, figsize=(6,4))
plt.subplots_adjust(wspace=0.4, hspace=0.6)
# (x1,y)
for irow in range(y.shape[0]):
    axs[0].scatter(x1[irow,:],y[irow,:])
axs[0].set_title('(x1,y)-scatter plot')
axs[0].set_xlabel('x1')
axs[0].set_ylabel('y')
# (x2,y)
for icol in range(y.shape[1]):
    axs[1].scatter(x2[:,icol],y[:,icol])
axs[1].set_title('(x2,y)-scatter plot')
axs[1].set_xlabel('x2')
axs[1].set_ylabel('y')
#
axs[0].grid(True)
axs[1].grid(True)
plt.show()<jupyter_output><empty_output><jupyter_text>三次元散布図の描画<jupyter_code># 三次元散布図の表示
fig = plt.figure()
ax = plt.axes(projection='3d')
ax.scatter(x1,x2,y)
ax.set_title('Demo')
ax.set_xlabel('x1')
ax.set_ylabel('x2')
ax.set_zlabel('y')
ax.set_xlim([minx1-0.1*rx1,maxx1+0.1*rx1])
ax.set_ylim([minx2-0.1*rx2,maxx2+0.1*rx2])
ax.view_init(elev=30,azim=30) # 見やすいように設定
ax.grid(True)
plt.show()<jupyter_output><empty_output><jupyter_text>線形重回帰の実行- sklearn.linear_model.LinearRegression を利用<jupyter_code>from sklearn.linear_model import LinearRegression

# 線形単回帰のインスタンス生成
regm = LinearRegression()

# 線形重回帰分析
Xdata = np.concatenate([x1.reshape(-1,1), x2.reshape(-1,1)],axis=1) # データ行列生成
ydata = y.reshape(-1,1) # 列ベクトル化
regm.fit(Xdata,ydata) # 最小２乗法
# 推定結果
beta0 = regm.intercept_[0]
beta1 = regm.coef_[0][0]
beta2 = regm.coef_[0][1]
print("β0 = {0:5.2f} ".format(beta0))
print("β1 = {0:5.2f} ".format(beta1))
print("β2 = {0:5.2f} ".format(beta2))<jupyter_output>β0 = 306.50 
β1 =  0.25 
β2 = -27140.07<jupyter_text>3次元散布図と曲面の表示<jupyter_code>nPoints = 10
fig = plt.figure()
ax = plt.axes(projection='3d')
xfit1,xfit2 = np.meshgrid(np.linspace(minx1,maxx1,nPoints),
                          np.linspace(minx2,maxx2,nPoints))
Xfit  = np.concatenate([xfit1.reshape(-1,1),xfit2.reshape(-1,1)],axis=1)
yfit  = regm.predict(Xfit).reshape(xfit1.shape)
ax.scatter(x1,x2,y)
ax.plot_wireframe(xfit1,xfit2,yfit, color = 'red')
ax.set_title('Demo')
ax.set_xlabel('x1')
ax.set_ylabel('x2')
ax.set_zlabel('y')
ax.set_xlim([minx1-0.1*rx1,maxx1+0.1*rx1])
ax.set_ylim([minx2-0.1*rx2,maxx2+0.1*rx2])
ax.view_init(elev=30,azim=30) # 見やすいように設定
ax.grid(True)
plt.show()<jupyter_output><empty_output><jupyter_text>基底展開重回帰モデルのパラメータの選択- RBF重回帰（K平均法を利用した2段階推定クラスを定義）- 交差検証正則化モデル（(sklearn.linear_model.LassoCV/RigdeCV）を利用<jupyter_code>"""GaussianFeaturesWithKmenasモジュール

    Copyright (c) 2020, <NAME>, All rights reserved.
"""
import numpy as np
import pandas as pd 
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression
from sklearn.pipeline import make_pipeline
from sklearn.base import BaseEstimator, TransformerMixin

class GaussianFeaturesWithKmeans(BaseEstimator, TransformerMixin):
    """K平均法を利用した2段階推定法によるガウス基底関数

        scikit-learn LinearRegression推定器に渡すためのガウス基底関数による
        特徴量への変換クラス    

    Attirbutes:
        nbfs (int): 基底関数の数
        centers_ (numpy.ndarray): 各ガウス関数の中心
        widths_ (numpy.ndarray): 各ガウス関数の幅

    References:

        * 小西貞則「多変量解析入門－線形から非線形へ－」岩波書店
        * Jake VanderPlas（菊池彰訳）「Pythonデータサイエンスハンドブック」オライリージャパン
    """
    def __init__(self, nbfs, width_factor=1.0, prekmeans=True):
        """GaussianFeaturesWithKMeansのコンストラクタ"""
        self.nbfs = nbfs
        self.width_factor = width_factor
        self.prekmeans = prekmeans

    @staticmethod
    def _gauss_basis(x, y, width, axis=None):
        arg = (x - y) / (width + 1e-300)
        return np.exp(-0.5 * np.sum(arg ** 2, axis))

    def fit(self, X, y=None):
        ndims = X.shape[1]
        if self.prekmeans: # K平均法を前処理に利用
            # 説明変数の標準化とK平均法の実行
            scaler = StandardScaler()         
            kmeans = KMeans(n_clusters=self.nbfs,random_state=0)
            X = scaler.fit_transform(X) 
            kmeans.fit(X) 
            # クラスタラベルの抽出
            labels = kmeans.predict(X).reshape(-1,1)            
            # 各クラスタの中心と幅を抽出
            clusters = pd.DataFrame(np.concatenate((labels,X),axis=1)).groupby([0])
            centers = scaler.inverse_transform(kmeans.cluster_centers_) 
            widths = (self.width_factor * scaler.scale_ * clusters.std(ddof=0)).to_numpy()
            #プライベート静的メソッド_gauss_basisでブロードキャスト計算するために成形
            if ndims == 1:
                self.centers_ = centers.reshape(-1,)
                self.widths_ = widths.reshape(-1,)
            else:
                self.centers_ = centers.reshape(-1,ndims,1).transpose(2,1,0) 
                self.widths_ = widths.reshape(-1,ndims,1).transpose(2,1,0)
        else: # 等間隔にガウス基底を配置（単変量のみ）
            if ndims == 1:
                centers = np.linspace(X.min(), X.max(), self.nbfs) 
                self.centers_ = centers
                self.widths_ = self.width_factor * (centers[1] - centers[0]) * np.ones(self.nbfs)
            else:
                raise Exception('prekmeansを True に設定してください．重回帰はK平均法による事前推定のみ対応しています．') 
        
        return self
        
    def transform(self, X):
        return self._gauss_basis(X[:,:, np.newaxis], self.centers_,
                                self.widths_, axis=1)<jupyter_output><empty_output><jupyter_text>正則化付き重回帰モデルの評価- リッジ交差検証(CV)モデル- ラッソ交差検証(CV)モデル<jupyter_code>from sklearn.linear_model import RidgeCV, LassoCV, Ridge, Lasso
from sklearn.preprocessing import PolynomialFeatures
from sklearn.metrics import mean_squared_error

# ラッソ回帰の選択
isLasso = True # False の場合はリッジ回帰を選択

# RFB重回帰の選択
isRBF = True # False の場合は多項式重回帰を選択
width_factor = 4.0 # ガウス関数の広さ

# 正則化法の選択
if isLasso: # ラッソ正則化交差検証(CV)モデル
    reg = LassoCV(alphas=np.logspace(-8,2,50))
else: # リッジ正則化交差検証(CV)モデル
    reg = RidgeCV(alphas=np.logspace(-8,2,50))

# 基底展開法の特徴量を選択
if isRBF: # RBF単回帰モデル
    # 評価する基底数の最大値（回帰分析が不安定になるときは値を減らして試す）
    maxdeg = 30
    phi = GaussianFeaturesWithKmeans(nbfs=maxdeg,width_factor=width_factor)
else: # 多項式単回帰モデル
    # 評価する基底数の最大値（回帰分析が不安定になるときは値を減らして試す）
    maxdeg = 10
    phi = PolynomialFeatures(degree=maxdeg,include_bias=False)
        
# 基底展開重回帰モデルのインスタンス生成
regmp = make_pipeline(phi,reg)
# 基底展開重回帰分析
regmp.fit(Xdata,ydata)


if isRBF: # RBF重回帰モデル
    nCoefs = maxdeg+1
else: # 多項式重回帰モデル
    nCoefs = int((maxdeg+1)*(maxdeg+2)/2) # (n+1)(n+2)/2

# 最適パラメータの推定結果
alpha = regmp[1].alpha_
print("λ = {0:5.2e} ".format(alpha))
w0 = regmp[1].intercept_
if not isLasso:
    w0 = w0[0]
print("w0 = {0:5.2f} ".format(w0))
for idx in range(0,nCoefs-1):
    if isLasso:
        w = regmp[1].coef_[idx]
    else:
        w = regmp[1].coef_[0][idx]
    print("w{0:d} = {1:5.2f} ".format(idx+1,w))<jupyter_output>λ = 1.39e-01 
w0 = 171.68 
w1 = 181.08 
w2 = -4.85 
w3 = -150.54 
w4 = -0.00 
w5 =  0.00 
w6 = 42.77 
w7 = -50.82 
w8 = 29.83 
w9 = -0.00 
w10 =  0.00 
w11 = -57.51 
w12 = 34.39 
w13 =  0.00 
w14 =  0.00 
w15 = -191.08 
w16 = 41.95 
w17 = 268.97 
w18 = -170.11 
w19 = -0.00 
w20 =  0.00 
w21 =  0.00 
w22 = 311.12 
w23 = -0.00 
w24 =  0.00 
w25 =  0.00 
w26 =  0.00 
w27 = 79.33 
w28 = -0.00 
w29 = 149.74 
w30 = -5.88<jupyter_text>3次元散布図と曲面の表示<jupyter_code># 回帰曲面をプロット
fig = plt.figure()
ax = plt.axes(projection='3d')
yfit = regmp.predict(Xfit).reshape(xfit1.shape)
ax.scatter(x1,x2,y)
ax.plot_wireframe(xfit1,xfit2,yfit, color = 'red')
ax.set_title("{0}... w/ {1}...".format(str(regmp[0])[:30],str(regmp[1])[:10]))
ax.set_xlabel('x1')
ax.set_ylabel('x2')
ax.set_zlabel('y')
ax.set_xlim([minx1-0.1*rx1,maxx1+0.1*rx1])
ax.set_ylim([minx2-0.1*rx2,maxx2+0.1*rx2])
ax.view_init(elev=30,azim=30) # 見やすいように設定
ax.grid(True)
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1519622307465329, 'fraction_numerical': 0.05586702075341792, 'mean_word_length': 4.0916374561842765, 'pattern_counts': {'":': 0, '<': 39, '<?xml version=': 0, '>': 49, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, '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': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '371640', 'license': 'MIT', 'n_tokens_mistral': 5378, 'n_tokens_neox': 5272, 'n_words': 634, 'path': 'example01_04.ipynb', 'repo_name': 'msiplab/EicEngLabIV', 'size': '10176'}