吴恩达作业ex5：Regularized Linear Regression and Bias v.s. Variance

可视化与加载数据

import numpy as np
import scipy.io as sio
from scipy.io import loadmat
import scipy.optimize as opt
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
1
2
3
4
5
6
7

d=loadmat('ex5data1.mat')
1

d
1

{'__header__': b'MATLAB 5.0 MAT-file, Platform: GLNXA64, Created on: Fri Nov  4 22:27:26 2011',
 '__version__': '1.0',
 '__globals__': [],
 'X': array([[-15.93675813],
        [-29.15297922],
        [ 36.18954863],
        [ 37.49218733],
        [-48.05882945],
        [ -8.94145794],
        [ 15.30779289],
        [-34.70626581],
        [  1.38915437],
        [-44.38375985],
        [  7.01350208],
        [ 22.76274892]]),
 'y': array([[ 2.13431051],
        [ 1.17325668],
        [34.35910918],
        [36.83795516],
        [ 2.80896507],
        [ 2.12107248],
        [14.71026831],
        [ 2.61418439],
        [ 3.74017167],
        [ 3.73169131],
        [ 7.62765885],
        [22.7524283 ]]),
 'Xtest': array([[-33.31800399],
        [-37.91216403],
        [-51.20693795],
        [ -6.13259585],
        [ 21.26118327],
        [-40.31952949],
        [-14.54153167],
        [ 32.55976024],
        [ 13.39343255],
        [ 44.20988595],
        [ -1.14267768],
        [-12.76686065],
        [ 34.05450539],
        [ 39.22350028],
        [  1.97449674],
        [ 29.6217551 ],
        [-23.66962971],
        [ -9.01180139],
        [-55.94057091],
        [-35.70859752],
        [  9.51020533]]),
 'ytest': array([[ 3.31688953],
        [ 5.39768952],
        [ 0.13042984],
        [ 6.1925982 ],
        [17.08848712],
        [ 0.79950805],
        [ 2.82479183],
        [28.62123334],
        [17.04639081],
        [55.38437334],
        [ 4.07936733],
        [ 8.27039793],
        [31.32355102],
        [39.15906103],
        [ 8.08727989],
        [24.11134389],
        [ 2.4773548 ],
        [ 6.56606472],
        [ 6.0380888 ],
        [ 4.69273956],
        [10.83004606]]),
 'Xval': array([[-16.74653578],
        [-14.57747075],
        [ 34.51575866],
        [-47.01007574],
        [ 36.97511905],
        [-40.68611002],
        [ -4.47201098],
        [ 26.53363489],
        [-42.7976831 ],
        [ 25.37409938],
        [-31.10955398],
        [ 27.31176864],
        [ -3.26386201],
        [ -1.81827649],
        [-40.7196624 ],
        [-50.01324365],
        [-17.41177155],
        [  3.5881937 ],
        [  7.08548026],
        [ 46.28236902],
        [ 14.61228909]]),
 'yval': array([[ 4.17020201e+00],
        [ 4.06726280e+00],
        [ 3.18730676e+01],
        [ 1.06236562e+01],
        [ 3.18360213e+01],
        [ 4.95936972e+00],
        [ 4.45159880e+00],
        [ 2.22763185e+01],
        [-4.38738274e-05],
        [ 2.05038016e+01],
        [ 3.85834476e+00],
        [ 1.93650529e+01],
        [ 4.88376281e+00],
        [ 1.10971588e+01],
        [ 7.46170827e+00],
        [ 1.47693464e+00],
        [ 2.71916388e+00],
        [ 1.09269007e+01],
        [ 8.34871235e+00],
        [ 5.27819280e+01],
        [ 1.33573396e+01]])}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111

数据包括：

• training set:X,y
• cross validation: Xval,yval
• test set:Xtest,ytest

X, y, Xval, yval, Xtest, ytest = map(np.ravel, [d['X'], d['y'], d['Xval'], d['yval'], d['Xtest'], d['ytest']])
1

df = pd.DataFrame({'water_level':X, 'flow':y})
df.head()
1
2

#可视化
sns.lmplot('water_level','flow',data=df,fit_reg=False)
plt.show()
1
2
3

E:\Anaconda\lib\site-packages\seaborn\_decorators.py:36: FutureWarning: Pass the following variables as keyword args: x, y. From version 0.12, the only valid positional argument will be `data`, and passing other arguments without an explicit keyword will result in an error or misinterpretation.
  warnings.warn(
1
2

Regularized linear regression cost function

X, Xval, Xtest = [np.insert(x.reshape(x.shape[0], 1), 0, np.ones(x.shape[0]), axis=1) for x in (X, Xval, Xtest)]
1

def cost(theta,X,y):
    #input:参数值theta,数据X，标签y
    #output：代价函数
    
    #样本个数
    m=X.shape[0]
    
    #计算代价函数
    inner=X@theta-y
    square_sum=inner.T@inner
    cost=square_sum/(2*m)
    
    return cost
1
2
3
4
5
6
7
8
9
10
11
12
13

theta = np.ones(X.shape[1])
cost(theta, X, y)
1
2

303.9515255535976
1

#正则化代价函数
def regularized_cost(theta,X,y,l=1):
    m=X.shape[0]
    
    regularized_term=(1/(2*m))*np.power(theta[1:],2).sum()
    return cost(theta,X,y)+regularized_term
1
2
3
4
5
6

Regularized linear regression gradient

def gradient(theta,X,y):
    #input:参数值theta，数据X，标签y
    #output：当前参数下梯度
    
    #样本个数
    m=X.shape[0]
    
    #计算代价函数
    grad=(X.T@(X@theta-y))/m
    
    return grad
1
2
3
4
5
6
7
8
9
10
11

gradient(theta, X, y)
1

array([-15.30301567, 598.16741084])
1

def regularized_gradient(theta, X, y, l=1):
    
    m=X.shape[0]
    
    regularized_term=theta.copy()
    regularized_term[0]=0#theta0不正则化
    
    regularized_term=(1/m)*regularized_term
    
    return gradient(theta, X, y) + regularized_term
1
2
3
4
5
6
7
8
9
10

regularized_gradient(theta, X, y)
1

array([-15.30301567, 598.25074417])
1

拟合数据

def linear_regression_np(X, y, l=1):

    # STEP1：初始化参数
    theta = np.ones(X.shape[1])
    
    # STEP2：调用优化算法拟合参数
    # your code here  (appro ~ 1 lines)
    res = opt.minimize(fun=regularized_cost,
                       x0=theta,
                       args=(X, y, l),
                       method='TNC',
                       jac=regularized_gradient,
                       options={'disp': True})
    return res

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

theta = np.ones(X.shape[0])

final_theta = linear_regression_np(X, y, l=0).get('x')
1
2
3

b=final_theta[0]
k=final_theta[1]

plt.scatter(X[:,1],y,label='training data',c='r')
plt.plot(X[:,1],X[:,1]*k+b,label='prediction')
plt.legend(loc=2)
plt.show()
1
2
3
4
5
6
7

bias-variance

training_cost, cv_cost = [], []
1

m = X.shape[0]
for i in range(1, m+1):
    # 计算当前样本的代价
    res = linear_regression_np(X[:i, :], y[:i], l=0)
    
    tc = regularized_cost(res.x, X[:i, :], y[:i], l=0)
    cv = regularized_cost(res.x, Xval, yval, l=0)
    
    # 把计算结果存储至预先定义的数组training_cost, cv_cost中
    training_cost.append(tc)
    cv_cost.append(cv)
1
2
3
4
5
6
7
8
9
10
11

plt.plot(np.arange(1,m+1),training_cost,label='training cost')
plt.plot(np.arange(1,m+1),cv_cost,label='cv cost')
plt.legend(loc=1)
plt.show()
1
2
3
4

多项式回归

#特征映射
def poly_features(x, power, as_ndarray=False):  #特征映射
    data = {'f{}'.format(i): np.power(x, i) for i in range(1, power + 1)}
    df = pd.DataFrame(data)

    return df.values if as_ndarray else df
1
2
3
4
5
6

def normalize_feature(df):
    """Applies function along input axis(default 0) of DataFrame."""
    return df.apply(lambda column: (column - column.mean()) / column.std())
1
2
3

def prepare_poly_data(*args, power):
    """
    args: keep feeding in X, Xval, or Xtest
        will return in the same order
    """
    def prepare(x):
        # 特征映射
        df = poly_features(x, power=power)

        # 归一化处理
        ndarr = normalize_feature(df).values

        # 添加偏置项
        return np.insert(ndarr, 0, np.ones(ndarr.shape[0]), axis=1)

    return [prepare(x) for x in args]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

X, y, Xval, yval, Xtest, ytest = map(np.ravel, [d['X'], d['y'], d['Xval'], d['yval'], d['Xtest'], d['ytest']])
1

poly_features(X, power=3)
1

画出线性回归

X_poly, Xval_poly, Xtest_poly= prepare_poly_data(X, Xval, Xtest, power=8)
1

def plot_learning_curve(X, y, Xval, yval, l=0):
# INPUT：训练数据集X,y，交叉验证集Xval，yval，正则化参数l
# OUTPUT：当前参数值下梯度
# TODO：根据参数和输入的数据计算梯度 
    
    # STEP1：初始化参数，获取样本个数，开始遍历
    training_cost, cv_cost = [], []
    m = X.shape[0]
    for i in range(1, m + 1):
        # STEP2：调用之前写好的拟合数据函数进行数据拟合
        # your code here  (appro ~ 1 lines)
        res = linear_regression_np(X[:i, :], y[:i], l=l)
        # STEP3：计算样本代价
        # your code here  (appro ~ 1 lines)
        tc = cost(res.x, X[:i, :], y[:i])
        cv = cost(res.x, Xval, yval)
        # STEP3：把计算结果存储至预先定义的数组training_cost, cv_cost中
        # your code here  (appro ~ 2 lines)
        training_cost.append(tc)
        cv_cost.append(cv)

        

    plt.plot(np.arange(1, m + 1), training_cost, label='training cost')
    plt.plot(np.arange(1, m + 1), cv_cost, label='cv cost')
    plt.legend(loc=1)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

plot_learning_curve(X_poly, y, Xval_poly, yval, l=0)
plt.show()
1
2

1