<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: python-dateutil>=2.1 in /usr/local/lib/python3.6/dist-packages (from matplotlib->japanize-matplotlib) (2.8.1)
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: six>=1.5 in /usr/local/lib/python3.6/dist-[...]<jupyter_text>Google Colaboratory でファイルを読み込む以下の方法があります。1. ファイルをアップロードする1. URL を指定する<jupyter_code># ファイルをアップロードする
#from google.colab import files
#uploaded_file = files.upload()
#filename = next(iter(uploaded_file))
#print(filename)

# 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>基底展開重回帰モデルのパラメータ数の選択- 多項式重回帰(sklearn.preprocessing.PolynomialFeaturessを利用)- RBF重回帰（K平均法を利用した2段階推定クラスを定義）<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>重回帰モデルの評価<jupyter_code>from sklearn.preprocessing import PolynomialFeatures
from sklearn.metrics import mean_squared_error

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

# 評価する基底数の最大値
maxdeg = 50

# モデル評価
deg = np.zeros(maxdeg,dtype=int)
aic = np.zeros(maxdeg)
regmp = [ [] for idx in range(maxdeg) ]
for idx in range(maxdeg):
    deg[idx] = idx+1

    # 基底展開法の特徴量を選択
    if isRBF: # RBF重回帰モデル
        phi = GaussianFeaturesWithKmeans(nbfs=deg[idx],width_factor=width_factor)
    else: #多項式重回帰モデル
        phi = PolynomialFeatures(degree=deg[idx],include_bias=False)
        
    # 基底展開重回帰モデルのインスタンス生成
    regmp[idx] = make_pipeline(phi,LinearRegression())
    # 基底展開重回帰分析
    regmp[idx].fit(Xdata,ydata)

    # AIC
    ypred = regmp[idx].predict(Xdata)
    mse = mean_squared_error(y.reshape(-1,1),ypred)
    if isRBF: # RBF重回帰モデル
        nCoefs = deg[idx]+1
    else: # 多項式重回帰モデル
        nCoefs = (deg[idx]+1)*(deg[idx]+2)/2 # (n+1)(n+2)/2
    aic[idx] = nSamples*np.log(2*np.pi*mse)+nSamples+2*(nCoefs+1)
    print("deg = {0:d}: AIC = {1:5.2e} ".format(deg[idx],aic[idx]))<jupyter_output>deg = 1: AIC = 2.24e+03 
deg = 2: AIC = 2.01e+03 
deg = 3: AIC = 1.75e+03 
deg = 4: AIC = 1.59e+03 
deg = 5: AIC = 1.59e+03 
deg = 6: AIC = 1.92e+03 
deg = 7: AIC = 2.22e+03 
deg = 8: AIC = 2.32e+03 
deg = 9: AIC = 2.39e+03 
deg = 10: AIC = 2.45e+03 
deg = 11: AIC = 2.50e+03 
deg = 12: AIC = 2.54e+03 
deg = 13: AIC = 2.59e+03 
deg = 14: AIC = 2.64e+03 
deg = 15: AIC = 2.70e+03 
deg = 16: AIC = 2.74e+03 
deg = 17: AIC = 2.76e+03 
deg = 18: AIC = 2.81e+03 
deg = 19: AIC = 2.85e+03 
deg = 20: AIC = 2.90e+03 
deg = 21: AIC = 2.95e+03 
deg = 22: AIC = 3.00e+03 
deg = 23: AIC = 3.04e+03 
deg = 24: AIC = 3.09e+03 
deg = 25: AIC = 3.15e+03 
deg = 26: AIC = 3.20e+03 
deg = 27: AIC = 3.26e+03 
deg = 28: AIC = 3.32e+03 
deg = 29: AIC = 3.38e+03 
deg = 30: AIC = 3.45e+03 
deg = 31: AIC = 3.51e+03 
deg = 32: AIC = 3.58e+03 
deg = 33: AIC = 3.65e+03 
deg = 34: AIC = 3.72e+03 
deg = 35: AIC = 3.81e+03 
deg = 36: AIC = 3.86e+03 
deg = 37: AIC = 3.94e+03 
deg = 38: AIC = 4.02e+03 
deg = 39: AIC = 4.10e[...]<jupyter_text>AICのプロット<jupyter_code>fig, ax = plt.subplots()
ax.plot(deg,aic)
plt.show()<jupyter_output><empty_output><jupyter_text>最小AICモデルの選択<jupyter_code># 最小のAICを与えるモデルのインデックス
optidx = np.argmin(aic)
# 最小のAICを与える次数
degopt = deg[optidx]
# 最適な基底展開重回帰のインスタンス抽出
regmpopt = regmp[optidx]
# 最適なパラメータ数
if isRBF: # RGB重回帰モデルのパラメータ数
    nCoefsOpt = int(degopt + 1)
else: # 多項式重回帰モデルのパラメータ数 
    nCoefsOpt = int((degopt+1)*(degopt+2)/2) # (n+1)(n+2)/2
# 最適パラメータの推定結果
w0 = regmpopt[1].intercept_[0]
print("w0 = {0:5.2f} ".format(w0))
for idx in range(0,nCoefsOpt-1):
    w = regmpopt[1].coef_[0][idx]
    print("w{0:d} = {1:5.2f} ".format(idx+1,w))<jupyter_output>w0 = -1.57 
w1 =  4.27 
w2 = -3.81 
w3 = -0.02 
w4 = -264.90 
w5 = -0.04 
w6 =  0.00 
w7 =  0.68 
w8 = -2.59 
w9 = -0.00 
w10 = -0.00 
w11 = -0.00 
w12 = 16.70 
w13 = -0.02 
w14 = -0.00 
w15 =  0.00 
w16 =  0.00 
w17 = -0.02 
w18 =  0.25 
w19 = -0.00 
w20 = -0.00<jupyter_text>3次元散布図と曲面の表示<jupyter_code>fig = plt.figure()
ax = plt.axes(projection='3d')
yfit = regmpopt.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>
{'fraction_non_alphanumeric': 0.15716051562775915, 'fraction_numerical': 0.06613102595797281, 'mean_word_length': 3.916232638888889, 'pattern_counts': {'":': 0, '<': 46, '<?xml version=': 0, '>': 57, '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': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '44674', 'license': 'MIT', 'n_tokens_mistral': 6101, 'n_tokens_neox': 5816, 'n_words': 775, 'path': 'example01_02.ipynb', 'repo_name': 'msiplab/EicEngLabIV', 'size': '11335'}