# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     name: python3
# ---

# + id="zU8oXGjshFjP" colab_type="code" colab={}
import numpy as np
# %matplotlib inline
import matplotlib.pyplot as plt
np.random.seed(0)
from statistics import mean

# + [markdown] id="x8buFpCZpHxN" colab_type="text"
# 今回はアルゴリズムの評価が中心の章なので,学習アルゴリズム実装は後に回し、sklearnを学習アルゴリズムとして使用する。

# + id="afwHFBQspYwV" colab_type="code" colab={}
import sklearn

# + [markdown] id="3DXxanUWrZLR" colab_type="text"
# 今回、学習に使うデータはsin関数に正規分布$N(\varepsilon|0,0.05)$ノイズ項を加えたデータを使う

# + id="w9qitzg_td9D" colab_type="code" colab={}
size = 100
max_degree = 11
x_data = np.random.rand(size) * np.pi * 2
var_data = np.random.normal(loc=0,scale=0.1,size=size)
sin_data = np.sin(x_data) + var_data

# + id="9JwIOkNav8lF" colab_type="code" outputId="b2f28b75-fdc4-407b-c4b5-4f81fa2e5af5" colab={"base_uri": "https://localhost:8080/", "height": 282}
plt.ylim(-1.2,1.2)
plt.scatter(x_data,sin_data)

# + [markdown] id="KJXHwwOgtfF4" colab_type="text"
#
# 学習用のアルゴリズムは多項式回帰を使います。

# + id="t4qyLuoQtecI" colab_type="code" colab={}
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
from sklearn.pipeline import Pipeline


# + [markdown] id="Pe_xlO-vDTv_" colab_type="text"
# 2.2.2:**MSE**:近似の良さの評価手法。
#
# $$MSE=\int (y(x;D) - h(x))^2p(x)dx=E\{(y(x;D)-h(x))^2\}$$

# + id="DXf_GKnqDS30" colab_type="code" colab={}
def MSE(y,t):
  return np.sum(np.square(y-t))/y.size


# + id="ZVYGA_lpR4PT" colab_type="code" outputId="2f41e964-3039-4250-c912-bcc1636e3d39" colab={"base_uri": "https://localhost:8080/", "height": 35}
MSE(np.array([10,3,3]),np.array([1,2,3]))


# + [markdown] id="WLknxmh7nUpV" colab_type="text"
# 2.2.1 (1)**ホールドアウト法**:
# 手元のデータを2つに分割し、片方をトレーニングに使い、片方をテストに使う手法。
# テストデータの数が必要

# + id="efYwVXAEoS47" colab_type="code" outputId="91cfd975-3c2d-466f-eac3-3ef06219894b" colab={"base_uri": "https://localhost:8080/", "height": 52}
# %%time
def holdout_method(x,y,per=0.8,value_func=MSE,degree=11):
  index = np.random.permutation(x.size)
  index_train,index_test = np.split(index,[int(x.size*per)])
  #plt.scatter(x_data[index_train],sin_data[index_train])
  test_score_list = []
  train_score_list = []
  for i in range(1,degree):
    pf = PolynomialFeatures(degree=i, include_bias=False)
    lr = LinearRegression()
    pl = Pipeline([("PF", pf), ("LR", lr)])
    pl.fit(x[index_train].reshape(-1,1), y[index_train])

    pred_y_test = pl.predict(x[index_test].reshape(-1,1))
    pred_y_train = pl.predict(x[index_train].reshape(-1,1))
    score_train = value_func(pred_y_train,y[index_train])
    score_test = value_func(pred_y_test,y[index_test])
    train_score_list.append(score_train)
    test_score_list.append(score_test)

  return train_score_list,test_score_list


# + id="vuWyhQp3LOqW" colab_type="code" outputId="14ce2f29-651c-4a27-d3d4-fdadbfbb6ff5" colab={"base_uri": "https://localhost:8080/", "height": 282}
hold_train_score_list,hold_test_score_list = holdout_method(x_data,sin_data,degree=max_degree)
plt.plot(np.array(range(1,max_degree)),np.array(hold_train_score_list),color='b')
plt.plot(np.array(range(1,max_degree)),np.array(hold_test_score_list),color='r')


# + [markdown] id="aCBHmr_o5Fpd" colab_type="text"
# (2)**交差確認法**:手元の各クラスをn分割して、n-1のグループで学習して、残りの１つのグループのデータでテストをし、その平均を誤り率とした性能評価を行う。

# + id="_iz_4m5f48ox" colab_type="code" colab={}
def cross_validation(x,y,value_func=MSE,split_num=5,degree=1):
  assert x.size % split_num==0,"You must use divisible number"
  n = x.size / split_num
  train_scores =[]
  test_scores =[]
  for i in range(split_num):
    indices = [int(i*n),int(i*n+n)]
    train_x_1,test_x,train_x_2=np.split(x,indices)
    train_y_1,test_y,train_y_2=np.split(y,indices)

    train_x = np.concatenate([train_x_1,train_x_2])
    train_y = np.concatenate([train_y_1,train_y_2])
    
    pf = PolynomialFeatures(degree=degree, include_bias=False)
    lr = LinearRegression()
    pl = Pipeline([("PF", pf), ("LR", lr)])
    pl.fit(train_x.reshape(-1,1), train_y)

    pred_y_test = pl.predict(np.array(test_x).reshape(-1,1))
    pred_y_train = pl.predict(np.array(train_x).reshape(-1,1))

    score_train = value_func(pred_y_train,train_y)
    #print(score_train)
    score_test = value_func(pred_y_test,test_y)
    #print(len(test_y))
    train_scores.append(score_train)
    test_scores.append(score_test)

  return mean(train_scores),mean(test_scores)


# + id="tYybB58UlhgR" colab_type="code" outputId="e89c7fea-6243-498f-d345-67a755af7772" colab={"base_uri": "https://localhost:8080/", "height": 282}
cross_test_score_list = []
cross_train_score_list = []
for i in range(1,max_degree):
  tra,tes = cross_validation(x_data,sin_data,degree=i)
  cross_train_score_list.append(tra)
  cross_test_score_list.append(tes)
plt.plot(np.array(range(1,max_degree)),np.array(cross_train_score_list),color='b')
plt.plot(np.array(range(1,max_degree)),np.array(cross_test_score_list),color='r')


# + [markdown] id="M3zr-OsK6vUk" colab_type="text"
# (3)**一つ抜き法**:交差確認法の特別な場合で、データ数＝グループの数としたものである。

# + id="r5oFd8dN5BWN" colab_type="code" colab={}
def leave_one_out(x,y,value_func=MSE,size=size,degree=1):
  return cross_validation(x,y,value_func,split_num=size,degree=degree)


# + id="P_NBcoykyOvL" colab_type="code" outputId="e932532e-4ff4-4720-dc0a-9c598eeb153f" colab={"base_uri": "https://localhost:8080/", "height": 282}
leave_test_score_list = []
leave_train_score_list = []
for i in range(1,max_degree):
  tra,tes = leave_one_out(x_data,sin_data,degree=i)
  leave_train_score_list.append(tra)
  leave_test_score_list.append(tes)
plt.plot(np.array(range(1,max_degree)),np.array(leave_train_score_list),color='b')
plt.plot(np.array(range(1,max_degree)),np.array(leave_test_score_list),color='r')

# + id="MfMN1SM_0cjh" colab_type="code" outputId="2f8d9b53-37a1-4ab0-dd2c-10e325413130" colab={"base_uri": "https://localhost:8080/", "height": 282}
plt.plot(np.array(range(1,max_degree)),np.array(hold_train_score_list),color='y')
plt.plot(np.array(range(1,max_degree)),np.array(hold_test_score_list),color='m')

plt.plot(np.array(range(1,max_degree)),np.array(cross_train_score_list),color='k')
plt.plot(np.array(range(1,max_degree)),np.array(cross_test_score_list),color='c')

plt.plot(np.array(range(1,max_degree)),np.array(leave_train_score_list),color='b')
plt.plot(np.array(range(1,max_degree)),np.array(leave_test_score_list),color='r')


# + [markdown] id="AcQKagJbCFnv" colab_type="text"
# (4)**ブートストラップ法**:N個の復元抽出をしてブートストラップサンプルを作り、そこから
#
# $bias=\varepsilon(N^*,N^*)-N(N^*,N)$
# を推定して、それをいくつか計算してその平均でバイアスを推定する。
# その推定値を$\overline{bias}$として、その推定値を
#
# $\varepsilon = \varepsilon(N,N)-\overline{bias}$
# とする。

# + id="dSu3P_fzCCVB" colab_type="code" colab={}
def bootstrap(x,y,value_func=MSE,trial=50,degree=1):
  biases=[]
  for i in range(trial):
    boot_ind = np.random.choice(range(x.size),size=x.size,replace=True)
    pf = PolynomialFeatures(degree=degree, include_bias=False)
    lr = LinearRegression()
    pl = Pipeline([("PF", pf), ("LR", lr)])
    pl.fit(x[boot_ind].reshape(-1,1), y[boot_ind])

    pred_y_boot = pl.predict(x[boot_ind].reshape(-1,1))
    pred_y_base = pl.predict(x.reshape(-1,1))

    score_boot = value_func(pred_y_boot,y[boot_ind])
    #print(score_train)
    score_base = value_func(pred_y_base,y)
    bias = score_base - score_boot
    #print(bias)
    biases.append(bias)

  pf = PolynomialFeatures(degree=degree, include_bias=False)
  lr = LinearRegression()
  pl = Pipeline([("PF", pf), ("LR", lr)])
  pl.fit(x.reshape(-1,1), y)

  pred_y_base = pl.predict(x.reshape(-1,1))
  score_base = value_func(pred_y_base,y)
  return score_base + mean(biases)


# + id="t2-ylH0gjjh4" colab_type="code" outputId="362a7bc4-c7aa-4adf-939c-3fedd8687def" colab={"base_uri": "https://localhost:8080/", "height": 282}
boot_score_list = []
for i in range(1,max_degree):
  boot_score = bootstrap(x_data,sin_data,degree=i)
  boot_score_list.append(boot_score)
plt.plot(np.array(range(1,max_degree)),np.array(boot_score_list),color='b')

StarcoderdataJupyterScriptsDedupFiltered