回帰分析においての損失関数の出し方(ridge)と交差検証の解釈について

Question

Pythonにて,Ridge回帰などを行い,その損失関数(MSE)の値を算出したいのですが,
その見方の方法がわからないためご教授していただけると幸いです.

また,Ridge回帰や線形単回帰,他項回帰で得られた,モデルの係数を
どのように参照したら良いかを教えていただけると幸いです.

また,付随する質問とはなりますが,
その損失関数をk-交差検証でモデルを評価した際に,
どのような判断をすれば良いかの質問を下部でしているため,お答えしていただけると幸いです.

何卒よろしくお願いいたします.

以下はコードと説明です.

# import library
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression
from sklearn.linear_model import Ridge
from sklearn.preprocessing import PolynomialFeatures

from sklearn.model_selection import cross_val_score
from sklearn.metrics import mean_squared_error
import statistics

# generate a true function
def func(x):
  return 100*x+10*x**3+0.001*x**15

# generate sample datas {(y,x)}_30
np.random.seed(0)
x = np.random.rand(30) * 2
noise = (np.random.rand(30) -0.5)* 100
true_y = func(x)
y = true_y + noise

# visualize true function and sample data
plt.scatter(x,y)
xx=np.arange(0,2,0.01)
yy = func(xx)
plt.plot(xx, yy)

x=x.reshape(-1,1)
# generate models

# degree=1
k_1=500
model1=LinearRegression()
model2=Ridge(alpha=k_1)

# degree = 20
k_20=500
model3=LinearRegression()
model4=Ridge(alpha=k_20)
base_x_20=PolynomialFeatures(degree=20,include_bias=False)
x_20=base_x_20.fit_transform(x)

model1.fit(x,y)
model2.fit(x,y)
model3.fit(x_20,y)
model4.fit(x_20,y)

# visualize model
# visualize model
def model_conf(model,x):
  y=[model.intercept_ for _ in range(len(x))]
  for j in range(len(x)):
    for i in range(len(model.coef_)):
      y[j]+=model.coef_[i]*((x[j][0])**(i+1))
  return y

plt.figure(figsize=(10, 7)) 

X_plt=np.arange(0,2,0.01).reshape(200,1)

y_pred = model1.predict(X_plt)
y_pred2 = model2.predict(X_plt)

line_20=base_x_20.transform(X_plt)
y_pred3=model3.predict(line_20)
y_pred4=model4.predict(line_20)

plt.scatter(x, y, color='blue', label='data')
plt.plot(X_plt, func(X_plt), color='purple', linestyle='-', label='True Function', lw=5)
plt.plot(X_plt, y_pred, color='red', linestyle=':', label='LinearRegression', lw=3)
plt.plot(X_plt, y_pred2, color='black', linestyle='--', label='Line Ridge(α={})'.format(k_1), lw=3)
#plt.plot(X_plt, model_conf(model3,X_plt), color='purple', linestyle="dashdot", label='Polynomial regression(20)', lw=5)
plt.plot(X_plt, y_pred3, color='green', linestyle="dashdot", label='Polynomial regression(20)', lw=3)
plt.plot(X_plt, y_pred4, color='pink', linestyle='-', label='Poly Ridge(α={})'.format(k_1), lw=3)
plt.legend()

plt.ylim([-50,300])
plt.show()

#損失関数を算出する.(修正 4/21, 17:29)
ans_lst=[]
for i in range(5):
  ans=0
  if i==0:
    #score=cross_val_score(model1,x,y,cv=10)
    ans=statistics.mean(list(map(lambda t: t**2, (y-model1.predict(x)).tolist())))
  elif i==1:
    model2.coef_
    ans=statistics.mean(list(map(lambda t: t**2, (y-model2.predict(x)).tolist())))
    ans+=sum(list(map(lambda t: t**2, model2.coef_)))
  elif i==2:
    ans=statistics.mean(list(map(lambda t: t**2, (y-model3.predict(x_20)).tolist())))
  elif i==3:
    ans=statistics.mean(list(map(lambda t: t**2, (y-model4.predict(x_20)).tolist())))
    ans+=sum(list(map(lambda t: t**2, model4.coef_)))
  else:
    #試しに真の分布のもののエラー関数を見てみる.
    ans=statistics.mean(list(map(lambda t: t**2, (y-func(x).reshape(1,-1)[0]).tolist())))
    ans_lst.append(ans)
    break
  ans_lst.append(ans)

"""
[878.1943851888309,
 6835.932288142904,
 165.22430413762072,
 1480.3156315725134,
 660.6537337438299]
"""

ans_lst=[878.1943851888309,
6829.205481962394,
165.22430413762072,
1472.468473432889,
660.6537337438299] <- ridge項を入れる前.

こちら,最後をみると明らかに過学習しているmodel3が
最も損失関数が低いため良いモデルとなってしまいました.

そのため,k-交差検証を用いて,それぞれのモデルを評価して,
その値が小さいモデルを判断します.

# generate models
models=[]

k_1=500
# degree=1
models.append(LinearRegression())
models.append(Ridge(alpha=k_1))

k_20=500
# degree = 20
models.append(LinearRegression())
models.append(Ridge(alpha=k_20))
base_x_20=PolynomialFeatures(degree=20,include_bias=False)
x_20=base_x_20.fit_transform(x)

# generate models
models=[]

# degree=1
k_1=500
models.append(LinearRegression())
models.append(Ridge(alpha=k_1))

# degree = 20
k_20=500
models.append(LinearRegression())
models.append(Ridge(alpha=k_20))
base_x_20=PolynomialFeatures(degree=20,include_bias=False)
x_20=base_x_20.fit_transform(x)

# evaluate model
scores=[]
for i in range(4):
  if i < 2 :
    scores.append(-np.mean(cross_val_score(models[i],x,y,scoring='neg_mean_squared_error',cv=10)))
  else:
    scores.append(-np.mean(cross_val_score(models[i],x_20,y,scoring='neg_mean_squared_error',cv=10)))
# 本当はridgeの部分はそれを含めたscoreを算出したいです.

"""
scores=[1050.7385117716492, 7591.142086205172, 5966152468.751239, 2145.3345481477745]
"""

この場合,models[0]を用いると良いとなりますが,結局最初のモデルを用いて
全てを訓練データとして得られたものを係数として答えを出せば良いでしょうか?
それとも,k-cross_validationで得られるそれぞれの係数を平均した方が良いのでしょうか?
こちらもご回答していただけると幸いです.

Answer 1

MSEの算出にはsklearn.metrics.mean_squared_errorを使います。真値と予測値を引数に与えればMSEを返してくれます。

係数は（存在すれば）coef_アトリビュートでアクセスできます。リファレンスのAttributes
の欄を確認してみてください。たとえばRidgeであれば、

coef_ndarray of shape (n_features,) or (n_targets, n_features)
Weight vector(s).

intercept_float or ndarray of shape (n_targets,)
Independent term in decision function. Set to 0.0 if fit_intercept = False.

sklearn.linear_model.Ridge — scikit-learn 0.24.1 documentation

といった項目が用意されています。

交差検証については、最良とされたモデルを全データで学習させるのが無難なやり方です。

---

なお、これを損失関数とは言いません。単に結果を評価しているだけなので、「MSEを評価指標に用いて評価した」等と書きます。

もしかしたら、ほんとうに最適化の過程で計算される損失関数（最小化される目的関数）の値を取得したいのかもしれませんが、scikit-learnでは残念ながら取得できません。

Comments

[unser]

ご回答いただきありがとうございます.
このような関数を用意すると確かに,coefで参照できていることが確認できました.
def model_conf(model,x):
y=[model.intercept_ for _ in range(len(x))]
for j in range(len(x)):
for i in range(len(model.coef_)):
y[j]+=model.coef_[i]*((x[j][0])**(i+1))
return y

確かに,損失関数というよりスコア関数の位置付けですね.ありがとうございます.

また,Ridge回帰において,k-crossvalidationの評価に関してscoringをneg_mean_squared_errorから何かに変更する必要があると思うのですが,そちらの方法はご存知でしょうか.

[hayataka2049]

＞また,Ridge回帰において,k-crossvalidationの評価に関してscoringをneg_mean_squared_errorから何かに変更する必要があると思うのですが,そちらの方法はご存知でしょうか.
特にその必要はないと思います。というか、他のモデルと違う指標を使ってしまうと、比較で困ります。

[unser]

なるほど.
つまり
目的関数は罰則項を含めた形で用いるけれど,
評価する際の関数に関しては
k-交差検証の場合は罰則項を含めずにやるということよろしいですか?

また今回Ridgeによる罰則をないモデルも用意しましたが,
今後はメインでRidge回帰同士で比較する場合を主に考慮したいためRidgeの罰則項を含めた評価実装も教えていただけると幸いです...

[hayataka2049]

目的関数と評価指標は独立に選べます。
（というか、「異なる目的関数を用いるモデル同士とか、そもそも目的関数という考え方を積極的に持たないようなモデルにおいて、優劣を比較したい」とかのために、評価指標を使う訳ですが）

Ridgeの罰則項を含めた評価をやる場合、
1.今回やっているのと同様に、罰則項のハイパーパラメータが異なる複数のRidgeモデルを作成する
2.GridSearchCVを使ってチューニングする
3.RidgeCVを使ってチューニングする

2が素直です。また、3はスマートかもしれません。
https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.GridSearchCV.html
https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.RidgeCV.html

[unser]

ありがとうございます.
確かに,独立に選べるべきですね!

確かに実用的にはGridSearchなどが良いですよね.
ベイズ最適化などでも実装したいと考えているので,その時わからなければまたこのサイトにて質問いたします!

ご丁寧にありがとうございます！

[jbpb0]

質問者さん

RidgeCVについて、このあたりが参考になるかもしれません
https://panda-clip.com/ridgecv/
http://www.science.smith.edu/~jcrouser/SDS293/labs/lab10-py.html