课程视频地址：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版作业4。

1. Neural Networks: MNIST image classfication

(a)(b)(c)前向传播部分很简单，反向传播部分的推导请参考CS231作业1，其中如下代码

#计算第二层梯度
p3 = np.exp(scores) / p2.reshape(-1, 1)
p3[np.arange(N), y] -= 1
dW2 = X1.T.dot(p3) / N + 2 * reg * W2
db2 = np.sum(p3, axis=0) / N
grads["W2"] = dW2
grads["b2"] = db2

对应于这里的

N2, D2 = y.shape
t2 = y - labels
db2 = np.sum(t2, axis=0) / N2
dW2 = h.T.dot(t2) / N2
if Lambda != 0:
	dW2 += 2 * Lambda * W2
dX2 = t2.dot(W2.T)

注意CS231的习题和这里最大的不同为前者的y是数值，即$0,1,…,9$的形式，而后者的则是one-hot的形式，所以以下两句代码作用相同：

#CS231
p3[np.arange(N), y] -= 1

#CS229
t2 = y - labels

在第一层中，只要知道sigmoid函数的导数为

即可，这部分的代码为：

#第一层梯度
N2, D2 = data.shape
t1 = h * (1 - h)
dW1 = data.T.dot(t1 * dX2) / N2
if Lambda != 0:
	dW1 += 2 * Lambda * W1
db1 = np.sum(t1, axis=0) / N2

完整的代码见my_nn_starter.py文件。

2. EM Convergence

首先回顾定义以及不等式：

如果$\theta$最终收敛到$\theta’$，那么

所以

3. PCA

首先计算$f_u(x)$，设

带入原式可得

关于$\alpha $求导，并令导数为$0$可得

所以

所以我们的目标为求解如下问题

记

那么

所以

注意前一项是常数，所以原问题可以化为

因此可以构造拉格朗日乘子：

求梯度可得

令上式为$0$，那么

这说明$u$是$X^TX$的特征值，并且

所以$u$是$X^TX$最大特征值对应的特征向量，即第一主成分。

4. Independent components analysis

报错参考：安装sounddevice报错的解决方案。

这里更新公式为：

其中

注意我们有

所以恢复数据的公式为

代码如下：

# -*- coding: utf-8 -*-
"""
Created on Fri Mar 29 13:49:57 2019

@author: qinzhen
"""

### Independent Components Analysis
###
### This program requires a working installation of:
###
### On Mac:
###     1. portaudio: On Mac: brew install portaudio
###     2. sounddevice: pip install sounddevice
###
### On windows:
###      pip install pyaudio sounddevice
###

import sounddevice as sd
import numpy as np

Fs = 11025

def normalize(dat):
    return 0.99 * dat / np.max(np.abs(dat))

def load_data():
    mix = np.loadtxt('mix.dat')
    return mix

def play(vec):
    sd.play(vec, Fs, blocking=True)

def sigmoid(x):
    return 1 / (1 + np.exp(-x))

def unmixer(X):
    M, N = X.shape
    W = np.eye(N)

    anneal = [0.1, 0.1, 0.1, 0.05, 0.05, 0.05, 0.02, 0.02, 0.01, 0.01,
              0.005, 0.005, 0.002, 0.002, 0.001, 0.001]
    print('Separating tracks ...')
    ######## Your code here ##########
    
    for alpha in anneal:
        #打乱数据
        np.random.permutation(X)
        for i in range(M):
            x = X[i, :].reshape(-1, 1)
            WX = W.dot(x)
            grad = (1 - 2 * sigmoid(WX)).dot(x.T) + np.linalg.inv(W.T)
            W += alpha * grad
    ###################################
    return W

def unmix(X, W):
    S = np.zeros(X.shape)

    ######### Your code here ##########
    S = X.dot(W.T)
    ##################################
    return S

X = normalize(load_data())

for i in range(X.shape[1]):
    print('Playing mixed track %d' % i)
    play(X[:, i])

W = unmixer(X)
S = normalize(unmix(X, W))

for i in range(S.shape[1]):
    print('Playing separated track %d' % i)
    play(S[:, i])

5. Markov decision processes

(a)$\forall \pi$，定义

那么

所以

因此

注意由定义我们有

由对称性，不妨设

如果记

从而

因此

类似的，如果

那么记

依然可得

所以无论那种情形，我们都有

对左边关于$s$取最大值可得

(b)如果存在$V$使得

那么任取同样满足上述条件的$V_0$，即

那么

因此

所以

6.Reinforcement Learning: The inverted pendulum

首先是初始化：

#价值函数
value_function = np.random.uniform(0, 0.1, NUM_STATES)
#计数
tran_cnt = np.zeros((NUM_STATES, NUM_STATES, NUM_ACTIONS))
#概率
tran_prob = np.ones((NUM_STATES, NUM_STATES, NUM_ACTIONS)) / NUM_STATES
#记录奖励
reward_value = np.zeros(NUM_STATES)
#记录奖励总数
reward_cnt = np.zeros(NUM_STATES)
#奖励
state_reward = np.zeros(NUM_STATES)

第一步初始化价值函数为均匀分布，第二步记录用行动$a$从状态$b$到状态$c$的次数，第三步初始化状态转移概率$P_{sa}$，第四步记录累计奖励，第五步统计奖励次数，第六步是初始化奖励函数$R(s)$。

然后是对行动作出选择：

#期望收益，比较其大小
s0 = np.sum(tran_prob[state, :, 0] * value_function)
s1 = np.sum(tran_prob[state, :, 1] * value_function)
if s0 > s1:
	action = 0
elif s0 < s1:
	action = 1
else:
	action = np.random.randint(0, 2)

上述代码第一步是计算

根据这一项选择采取哪种行动。

其次是每次行动后对统计量进行更新：

tran_cnt[state, new_state, action] += 1
reward_value[new_state] += R
reward_cnt[new_state] += 1

接着在到达最终状态后，更新价值函数，状态转移概率：

#统计在某个状态采取某个动作的次数
for i in range(NUM_ACTIONS):
	for j in range(NUM_STATES):
		#在状态j采取行动i的总数
		num = np.sum(tran_cnt[j, :, i])
		if num > 0:
			 tran_prob[j, :, i] = tran_cnt[j, :, i] / num
#更新奖励
state_reward[reward_cnt>0] = reward_value[reward_cnt>0] / reward_cnt[reward_cnt>0]

最后一步利用值迭代更新：

  • 对每个状态$s$，初始化$V(s):=0$

  • 重复直到收敛{

    • 对每个状态，更新$V(s):=R(s) + \max_{a\in A}\gamma \sum_{s’\in S}
      P_{sa}(s’) V(s’)$

    }

代码如下

iterations = 0
while True:
	#值迭代更新后的值
	v0 = GAMMA * tran_prob[:, :, 0].dot(value_function) + state_reward
	v1 = GAMMA * tran_prob[:, :, 1].dot(value_function) + state_reward
	v = np.c_[v0, v1]
	#取最大值
	value_function_new = np.max(v, axis=1)
	#计算更新幅度
	delta = np.linalg.norm(value_function_new - value_function)
	#更新价值函数
	value_function = np.copy(value_function_new)
	iterations += 1
	#更新幅度变小，则停止循环
	if delta < TOLERANCE:
		break

#更新一次收敛的计数
if iterations == 1:
	consecutive_no_learning_trials += 1
else:
	consecutive_no_learning_trials = 0