CS229 2017版作业1

课程视频地址:http://open.163.com/special/opencourse/machinelearning.html

课程主页:http://cs229.stanford.edu/

更具体的资料链接:https://www.jianshu.com/p/0a6ef31ff77a

作业地址:https://github.com/Doraemonzzz/CS229

参考资料:https://github.com/zyxue/stanford-cs229

这部分回顾2017版作业1。

1.Logistic regression

(a)首先回顾$J(\theta)$的定义

注意$g^{‘}(z) = g(z)(1-g(z))$,利用这点来求$\frac{\partial h_θ(x) }{\partial \theta_k}$

利用$\frac{\partial h_θ(x) }{\partial \theta_k}$来求$\frac{\partial J(\theta)}{\partial \theta_k}$

这个形式方便我们求二阶偏导,继续利用$\frac{\partial h_θ(x) }{\partial \theta_k}$来求$\frac{\partial^2 J(\theta)}{\partial \theta_l\partial \theta_k}$

接着作以下记号

所以Hessian矩阵$H$可以表达为如下形式

为了(b)题需要,这里也将$\nabla J(\theta)$表示出来

现在任取$z\in \mathbb R^n$,记$t=X^Tz$,那么

注意$h_θ(y^{(i)}x^{(i)}) \in [0,1]$,所以$h_θ(y^{(i)}x^{(i)})(1-h_θ(y^{(i)}x^{(i)}) \ge 0$,从而

从而$H$为半正定矩阵。

(b)(c)

这里解释下计算步骤,第一步是读取数据并增加一个截距项,由以下两个函数完成。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
%matplotlib inline
import numpy as np
#import matplotlib as mpl
#mpl.use('Agg')
import matplotlib.pyplot as plt
from numpy.linalg import inv
#from __future__ import division


def load_data():
X = np.genfromtxt('logistic_x.txt')
Y = np.genfromtxt('logistic_y.txt')
return X, Y

#增加截距项
def add_intercept(X_):
m, n = X_.shape
X = np.zeros((m, n + 1))
################
ones = np.ones((m, 1))
X = np.append(ones, X_, axis = 1)
################
return X

第二步是利用刚刚的公式计算梯度以及Hessian矩阵。

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
#利用之前所述的公式计算
def calc_grad(X, Y, theta):
m, n = X.shape
grad = np.zeros(theta.shape)

##############
Y_ = Y.reshape([-1, 1])
d1 = (X*Y_).dot(theta)
h = 1 / (1 + np.exp(-d1))
S = (1 - h) * Y
grad = -1/m * (X.T).dot(S)

##############

return grad

def calc_hessian(X, Y, theta):
m, n = X.shape
H = np.zeros((n, n))

##############
Y = Y.reshape([-1, 1])
d1 = (X*Y).dot(theta)
h = 1 / (1 + np.exp(-d1))
S = np.diag(h * (1-h))
H = X.T.dot(S).dot(X)

#############

return H

这里还有两个辅助的函数,分别是计算$J(\theta)$和作图。

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
def calc_loss(X, Y, theta):
m, n = X.shape
loss = 0.

###########
Y = Y.reshape([-1, 1])
d1 = (X*Y).dot(theta)
h = 1 / (1 + np.exp(-d1))
loss = -1/m * np.sum(np.log(h))

###########

return loss

def plot(X, Y, theta):
#plt.figure()

############
x1 = X[Y>0][:, 1]
y1 = X[Y>0][:, 2]
x2 = X[Y<0][:, 1]
y2 = X[Y<0][:, 2]
#计算系数
theta = logistic_regression(X, Y)
Min = np.min(X[:, 1])
Max = np.max(X[:, 1])
x = np.array([Min, Max])
y = -(theta[0] + theta[1]*x)/theta[2]
plt.scatter(x1, y1)
plt.scatter(x2, y2)
plt.plot(x, y)
plt.title('Newton’s method for Logistic regression')
############

plt.savefig('ps1q1c.png')
return

最重要的一步就是利用如下公式迭代计算:

1
2
3
4
5
6
7
8
9
10
11
def logistic_regression(X, Y):
m, n = X.shape
theta = np.zeros(n)

############
H = calc_hessian(X, Y, theta)
grad = calc_grad(X, Y, theta)
theta -= inv(H).dot(grad)
############

return theta
1
2
3
4
X_, Y = load_data()
X = add_intercept(X_)
theta = logistic_regression(X, Y)
plot(X, Y, theta)

全部代码可以查看my_logistic.py这个文件。

2.Poisson regression and the exponential family

(a)

对比指数族的形式

可得

所以泊松分布为指数族。

(b)我们来计算$g(η) = \mathbb E[T (y); η]$

(c)根据GLM的性质可知$\eta = \theta^T x$,那么$\lambda=e^{\eta} =e^{\theta^T x}$

计算对数似然函数

关于$\theta_j$求偏导可得

此处求最大值,用随机梯度上升法,更新规则为

(d)$T(y)=y$,所以

因为$\eta = \theta^Tx$,所以

接下来只要证明$\frac{\partial a(η)}{\partial η}= h(x) = \mathbb E[y;\eta]$即可,利用$p(y; η) $为概率密度函数

两边关于$\eta$求偏导可得

所以

从而梯度上升法的更新规则为

3.Gaussian discriminant analysis

(a)先计算$P(x)$

利用贝叶斯公式计算$P(y|x)$,分$y=1,y=-1$计算

所以

可以看到,指数部分都有一样的式子,现在计算这个式子

所以

综合两式,$P(y|x) $可以写为

从而

(b)(c)

(c)是(b)的一般情形,所以直接处理(c)

先观察$P(x|y)$的形式,可以得到如下公式

接着计算$\text{log}P(x,y)$

对数似然函数为

关于$\phi$求梯度

关于$\mu_1,\mu_{-1}$求梯度

求$\Sigma$的时候利用一些技巧性,我们不求的$\Sigma$极大似然估计,而是求$\Sigma^{-1}$的极大似然估计,然后再求出$\Sigma $的极大似然估计,利用如下两个式子

那么

所以结论成立。

4.Linear invariance of optimization algorithms

(a)我们计算$∇_zg(z),∇_z^2g(z)$。

先计算$∇_zg(z),$

所以

接着计算$∇_z^2g(z)$

从而

接着利用数学归纳法来证明结论。

$n=0$时,

所以$n=0$时结论成立。假设$n\le i$时,$z^{(n)} = A^{-1}x^{(n)}$,那么$n=i+1$时

其中倒数第二步是因为$x^{(i)} = Az^{(i)}$,所以$n=i+1$时结论成立,从而牛顿法满足invariant to linear reparameterizations

(b)对于梯度下降法, 继续利用

假设$z^{(i)} = A^{-1} x^{(i)}$,那么

但是

所以$z^{(i+1)} $与$A^{-1}x^{(i+1)}$不相等,从而梯度下降法不满足invariant to linear reparameterizations

5.Regression for denoising quasar spectra

(a)

(i)

那么

(ii)

(iii)

似然函数为

对于固定的$\sigma^{(i)}$,最大化这个概率似然函数等价于最小化

记$w^{(i)}=\frac 1 { {(σ^{(i)})^2}}$,这个式子可以转化为

(b)

(i)

第一步还是读取数据以及增加截距。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
%matplotlib inline
import numpy as np
#import matplotlib as mpl
#mpl.use('Agg')
import matplotlib.pyplot as plt
from numpy.linalg import inv

def load_data():
train = np.genfromtxt('quasar_train.csv', skip_header=True, delimiter=',')
test = np.genfromtxt('quasar_test.csv', skip_header=True, delimiter=',')
wavelengths = np.genfromtxt('quasar_train.csv', skip_header=False, delimiter=',')[0]
return train, test, wavelengths

def add_intercept(X_):
X = None
#####################
X_ = X_.reshape((-1, 1))
m, n = X_.shape
ones = np.ones((m, 1))
X = np.append(ones, X_, axis = 1)

###################
return X

这一部分直接利用最小二乘法计算结果。

1
2
3
4
5
6
7
8
9
10
11
12
#利用最小二乘公式
def LR_smooth(Y, X_):
X = add_intercept(X_)
Y = Y.reshape((-1, 1))
yhat = np.zeros(Y.shape)
theta = np.zeros(2)
#####################
theta = inv(X.T.dot(X)).dot(X.T).dot(Y)
yhat = X.dot(theta)

#####################
return yhat, theta

作图函数

1
2
3
4
5
6
7
8
9
10
11
def plot_b(X, raw_Y, Ys, desc, filename):
plt.figure()
############
for i in range(len(Ys)):
plt.scatter(X, raw_Y, c='red', s=2, label='raw_data')
plt.plot(X, Ys[i], c='blue', label='Regression line')
plt.title(desc[i])
plt.show()

############
plt.savefig(filename)

产生结果

1
2
3
4
5
6
raw_train, raw_test, wavelengths = load_data()

## Part b.i
lr_est, theta = LR_smooth(raw_train[0], wavelengths)
print('Part b.i) Theta=[%.4f, %.4f]' % (theta[0], theta[1]))
plot_b(wavelengths, raw_train[0], [lr_est], ['Regression line'], 'ps1q5b1.png')
Part b.i) Theta=[2.5134, -0.0010]

png

(ii)利用a中的公式$\theta =( X^T WX )^{-1} X^TWy$计算

1
2
3
4
5
6
7
8
9
10
11
12
13
def LWR_smooth(spectrum, wavelengths, tau):
smooth_spectrum = np.array([])
###############
X = add_intercept(wavelengths)
Y = spectrum.reshape((-1, 1))
for i in range(len(wavelengths)):
w = np.exp(-(wavelengths - wavelengths[i])**2 / (2*tau**2))
W = np.diag(w)
theta = inv(X.T.dot(W).dot(X)).dot(X.T).dot(W).dot(Y)
yhat = theta.T.dot(X[i])
smooth_spectrum = np.append(smooth_spectrum, yhat)
###############
return smooth_spectrum

计算结果

1
2
3
## Part b.ii
lwr_est_5 = LWR_smooth(raw_train[0], wavelengths, 5)
plot_b(wavelengths, raw_train[0], [lwr_est_5], ['tau = 5'], 'ps1q5b2.png')

png

对不同的参数作图。

1
2
3
4
5
6
7
8
9
### Part b.iii
lwr_est_1 = LWR_smooth(raw_train[0], wavelengths, 1)
lwr_est_10 = LWR_smooth(raw_train[0], wavelengths, 10)
lwr_est_100 = LWR_smooth(raw_train[0], wavelengths, 100)
lwr_est_1000 = LWR_smooth(raw_train[0], wavelengths, 1000)
plot_b(wavelengths, raw_train[0],
[lwr_est_1, lwr_est_5, lwr_est_10, lwr_est_100, lwr_est_1000],
['tau = 1', 'tau = 5', 'tau = 10', 'tau = 100', 'tau = 1000'],
'ps1q5b3.png')

png

png

png

png

(c)

(i)

第一步是利用局部加权回归来使得数据更平滑一些。

1
2
3
4
5
6
7
8
9
10
def smooth_data(raw, wavelengths, tau):
smooth = None
################
smooth = []
for spectrum in raw:
smooth.append(LWR_smooth(spectrum, wavelengths, tau))
################
return np.array(smooth)

smooth_train, smooth_test = [smooth_data(raw, wavelengths, 5) for raw in [raw_train, raw_test]]

(ii)
这部分比较难懂,详细解释下各个步骤,第一步将数据分为左边和右边,波长小于1200的为left,大于1300的为right,利用如下函数完成这部分工作。

1
2
3
4
5
6
7
8
9
def split(full, wavelengths):
left, right = None, None
###############
indexl = np.argwhere(wavelengths == 1200)[0][0]
indexr = np.argwhere(wavelengths == 1300)[0][0]
left = full[:, :indexl]
right = full[:, indexr:]
###############
return left, right

第二步要计算距离,定义如下函数

1
2
3
4
5
6
def dist(a, b):
dist = 0
################
dist = ((a - b)**2).sum()
################
return dist
1
2
3
#### Part c.ii
left_train, right_train = split(smooth_train, wavelengths)
left_test, right_test = split(smooth_test, wavelengths)

最后一步是最复杂的,首先计算距离,这里的距离定义为

对于此题来说,$f_1(\lambda),f_2(\lambda)$分别对应了测试数据以及训练数据,利用下式计算距离矩阵。

1
2
3
4
d = (right_train - right_test)**2

#求和
d1 = d.sum(axis=1)

然后要找到距离最近的$k$个点,这里为$3$个点,我的思路是先对数据进行排序,然后找到前三个数据的索引

1
2
3
4
5
#找到排名前3的作为neighb_k(f_right)
tempd = d1.copy()
tempd.sort()
#找到索引
index = (d1==tempd[0])|(d1==tempd[1])|(d1==tempd[2])

接着计算$h$,为距离的最大值,根据题目中的公式对数据进行转换

h = d1.max()
d1 = d1/h

最后一步根据kernel以及题目中的公式计算即可

1
2
3
4
5
6
d1 = 1 - d1
#求lefthat
a = (d1[index].dot(left_train[index]))
b = d1[index].sum()

lefthat = a/b

以上全部综合起来就是如下函数:

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
def func_reg(left_train, right_train, right_test):
m, n = left_train.shape
#m=200,n=50
lefthat = np.zeros(n)
###########################
#right_train 200*300
#right_test 1*300
#left_train 200*50
#求题目中的d(f1,f2),先求每个点的距离,200*300矩阵
d = (right_train - right_test)**2
#按照行求和200*1
d1 = d.sum(axis=1)
#找到排名前3的作为neighb_k(f_right)
tempd = d1.copy()
tempd.sort()
#找到索引
index = (d1==tempd[0])|(d1==tempd[1])|(d1==tempd[2])
#h为d1中的最大值
h = d1.max()
d1 = d1/h
#ker (1-t)
d1 = 1 - d1
#求lefthat
a = (d1[index].dot(left_train[index]))
b = d1[index].sum()

lefthat = a/b
###########################
return lefthat

作图函数

1
2
3
4
5
6
7
8
9
10
#将左边右边区分开来,左边<1200,右边>=1300
def plot_c(Yhat, Y, X, filename):
plt.figure()
############
plt.plot(X[:50],Yhat)
plt.plot(X,Y)
plt.show()
#############
plt.savefig(filename)
return

产生结果

1
2
3
4
5
6
#### Part c.ii
left_train, right_train = split(smooth_train, wavelengths)
left_test, right_test = split(smooth_test, wavelengths)

train_errors = [dist(left, func_reg(left_train, right_train, right)) for left, right in zip(left_train, right_train)]
print('Part c.ii) Training error: %.4f' % np.mean(train_errors))
Part c.ii) Training error: 1.0664

(iii)

1
2
3
4
5
6
7
8
### Part c.iii
test_errors = [dist(left, func_reg(left_train, right_train, right)) for left, right in zip(left_test, right_test)]
print('Part c.iii) Test error: %.4f' % np.mean(test_errors))

left_1 = func_reg(left_train, right_train, right_test[0])
plot_c(left_1, smooth_test[0], wavelengths, 'ps1q5c3_1.png')
left_6 = func_reg(left_train, right_train, right_test[5])
plot_c(left_6, smooth_test[5], wavelengths, 'ps1q5c3_6.png')
Part c.iii) Test error: 2.7100

png

png

全部代码可以查看my_quasars.py这个文件。

本文标题:CS229 2017版作业1

文章作者:Doraemonzzz

发布时间:2019年03月01日 - 17:54:00

最后更新:2019年03月01日 - 18:02:32

原始链接:http://doraemonzzz.com/2019/03/01/CS229 2017版作业1/

许可协议: 署名-非商业性使用-禁止演绎 4.0 国际 转载请保留原文链接及作者。