简易机器学习代码（LR，Kmeans，NN，RNN）

Logistic Regression

特别需要注意的是 exp 和 log 的使用。

sigmoid 原始表达式为 1 / (1+exp(-z))，但如果直接使用 z=-710，会显示 overflow。因此对于 z<0 的情况，采用 exp(z) / (1 + exp(z)) ，这样一来，exp(-710) 就没问题了。这就是 scipy 包里的 expit 函数

log_logistic = log(sigmoid)，注意和 expit 函数是一致的，分情况讨论。

import numpy as np
from scipy.special import expit
from sklearn.utils.extmath import log_logistic

def predict(theta, x):
    return expit(x.dot(theta))

def compute_loss(y, yz):
    return - np.sum(log_logistic(yz))

def gradientdescent(x, y, theta, iterations=2000, lr=0.01):
    loss_list = []
    for i in range(iterations):
        yhat = predict(theta, x)
        delta = x.T.dot(yhat - y) / m
        loss = compute_loss(y, y * X.dot(theta))
        loss_list.append(loss)
        theta = theta - lr * delta
    return theta, loss_list

theta, loss_list = gradientdescent(X, y, np.zeros((n, 1)))

Kmeans

Kmeans 的本质就是 EM 算法，只不过是硬间隔而不是软间隔。首先初始化 K 个中心点，在 E 步，将样本分配到最近的中心点，在 M 步，选取新的中心点以最小化组内距离。

import numpy as np


def calc_dist(x1, x2):
    return sum([(x1[i] - x2[i])**2 for i in range(len(x1))])


# Assign samples to given centers
def E_step(X, cents):
    cent_dict = dict(zip(cents, [[] for _ in range(len(cents))]))
    for row in X:
        min_dist, best_cent = 1e10, None
        for cent in cent_dict:
            dist = calc_dist(row, cent)
            if dist < min_dist:
                min_dist = dist
                best_cent = cent
        cent_dict[best_cent] += [row.tolist()]
    return cent_dict
    

# Compute new centers
def M_step(cent_dict):
    new_cents = []
    for cent in cent_dict:
        new_cent = np.mean(np.array(cent_dict[cent]), axis=0)
        new_cents.append(tuple(new_cent))
    return new_cents


def Kmeans(X, K=3, max_iter=10):
    np.random.seed(1)
    inds = np.random.choice(len(X), K)
    init_cents = [tuple(X[i]) for i in inds]
    cents = init_cents
    for k in range(max_iter):
        cent_dict = E_step(X, cents)
        new_cents = M_step(cent_dict)
        move = sum([calc_dist(c1, c2) for c1, c2 in zip(cents, new_cents)])
        if move < 0.1:
            print('Converged in %s steps' % k)
            break
        cents = new_cents
    return cent_dict

Neural Network

注意 softmax 的计算，需要考虑到 exp 的 overflow。因此通常会在 softmax 分子分母同时乘上一个常数 C，log(C) = -max(z)，这就是 shift_score。

这里使用了 scipy 包里的 logsumexp，理由同 LR，logsumexp = log(sum(exp(z)))。

from scipy.special import logsumexp
import numpy as np



class Neural_Network:
    

    def __init__(self, n, h, c, std=1e-4):
        W1 = np.random.randn(n, h) * std
        b1 = np.zeros(h)
        W2 = np.random.randn(h, c) * std
        b2 = np.zeros(c)
        self.params = {'W1': W1, 'b1': b1, 'W2': W2, 'b2': b2}
        
        
    def forward_backward_prop(self, X, y):
        W1, b1 = self.params['W1'], self.params['b1']
        W2, b2 = self.params['W2'], self.params['b2']
        
        # forward prop
        hidden = X.dot(W1) + b1
        relu = np.maximum(0, hidden)
        scores = relu.dot(W2) + b2
        shift_scores = scores - np.max(scores, axis=1, keepdims=True)
        softmax = np.exp(shift_scores) / np.sum(np.exp(shift_scores), axis=1, keepdims=True)
        loss = - np.sum(y * (shift_scores - logsumexp(shift_scores, axis=1, keepdims=True))) / X.shape[0]
        
        # backward prop
        dscores = (softmax - y) / X.shape[0]
        drelu = dscores.dot(W2.T)
        dW2 = relu.T.dot(dscores)
        db2 = np.sum(dscores, axis=0)
        dhidden = (hidden > 0) * drelu
        dW1 = X.T.dot(dhidden)
        db1 = np.sum(dhidden, axis=0)
        
        grads = {'dW1': dW1, 'db1': db1, 'dW2': dW2, 'db2': db2}
        
        return loss, grads
        
    
    def train(self, X, y, lr=0.01, decay=0.95, iters=5000):
        loss_list, acc_list = [], []
        for it in range(iters):
            loss, grads = self.forward_backward_prop(X, y)
            loss_list.append(loss)
            self.params['W1'] -= lr * grads['dW1']
            self.params['b1'] -= lr * grads['db1']
            self.params['W2'] -= lr * grads['dW2']
            self.params['b2'] -= lr * grads['db2']
            
            if it % 100 == 0:
                yhat = self.predict(X)
                acc = np.sum(np.argmax(y, axis=1) == yhat) / X.shape[0]
                acc_list.append(acc)
                lr *= decay
                
        return loss_list, acc_list
    
    
    def predict(self, X):
        hidden = X.dot(self.params['W1']) + self.params['b1']
        relu = np.maximum(0, hidden)
        scores = relu.dot(self.params['W2']) + self.params['b2']
        yhat = np.argmax(scores, axis=1)
        return yhat

Recurrent Neural Network

import numpy as np



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

def softmax(x):
    ex = np.exp(x - np.max(x))
    return ex / ex.sum(axis=0)



class RNN:
    

    def __init__(self, na, nx, ny, m, seed=1):
        np.random.seed(seed)
        Waa = np.random.randn(na, na)
        Wax = np.random.randn(na, nx)
        Wya = np.random.randn(ny, na)
        ba = np.random.randn(na, 1)
        by = np.random.randn(ny, 1)
        self.a0 = np.random.randn(na, m)
        self.params = {'Waa': Waa, 'Wax': Wax, 'Wya': Wya, 'ba': ba, 'by': by}
    
    
    def RNN_cell_forward(self, xt, a_prev):
        """
        Inputs:
        xt -- Current input data, of shape (nx, m).
        a_prev -- Previous hidden state, of shape (na, m)
        
        Outputs:
        at -- Current hidden state, of shape (na, m)
        yt -- Current prediction, of shape (ny, m)
        """
        Waa, Wax, ba = self.params['Waa'], self.params['Wax'], self.params['ba']
        Wya, by = self.params['Wya'], self.params['by']
        
        at = tanh(Waa.dot(a_prev) + Wax.dot(xt) + ba)
        score = Wya.dot(at) + by
        yt = softmax(score)
        return at, yt
    
    
    def RNN_forward(self, X, y):
        """
        Inputs:
        X -- Input data for every time step, of shape (nx, m, Tx)
        y -- Target for every time step, of shape (ny, m, Tx)
        
        Outputs:
        a -- Hidden states for every time-step, of shape (n_a, m, T_x)
        yhat -- Predictions for every time-step, of shape (n_y, m, T_x)
        """
        a_prev = self.a0
        na, m = a_prev.shape
        ny = y.shape[0]
        Tx = X.shape[2]
        
        a = np.zeros((na, m, Tx))
        yhat = np.zeros((ny, m, Tx))
        loss = 0
        for t in range(Tx):
            a_next, yt = self.RNN_cell_forward(X[:, :, t], a_prev)
            yhat[:, :, t] = yt
            a[:, :, t] = a_next
            loss -= np.sum(np.log(yt.T.dot(y[:, :, t])))
            a_prev = a_next
            
        cache = (a, yhat)
        return loss, cache
    
    
    def RNN_cell_backward(self, dz, grads, cache):
        """
        Inputs:
        dz -- Gradient of loss with respect to score
        grads -- Dictionary contains all gradients
        cache -- Tuple contains xt, a_next, a_prev
        
        Outputs:
        grads -- Dictionary contains all gradients
        """
        xt, a_next, a_prev = cache
        Waa, Wax, ba = self.params['Waa'], self.params['Wax'], self.params['ba']
        Wya, by = self.params['Wya'], self.params['by']
        
        grads['dWya'] += dz.dot(a_next.T)
        grads['dby'] += np.sum(dz, axis=1, keepdims=True)
        da_y = Wya.T.dot(dz) 
        da_a = grads['da_prev']
        da_next = da_y + da_a      # da is computed based on two paths, from da_y and da_a.
        dtanh = (1 - a_next**2) * da_next
        grads['dWaa'] += dtanh.dot(a_prev.T)
        grads['da_prev'] = Waa.T.dot(dtanh)
        grads['dWax'] += dtanh.dot(xt.T)
        grads['dba'] += np.sum(dtanh, axis=1, keepdims=True)
        
        return grads
    
    
    def RNN_backward(self, X, y, cache):
        """
        Inputs:
        X -- Input data for every time step, of shape (nx, m, Tx)
        y -- Target for every time step, of shape (ny, m, Tx)
        cache -- Tuple from RNN_forward, contains a, yhat
        
        Outputs:
        grads -- Dictionary contains all gradients
        a -- Hidden states for every time-step, of shape (n_a, m, T_x)
        """
        a, yhat = cache
        Waa, Wax, ba = self.params['Waa'], self.params['Wax'], self.params['ba']
        Wya, by = self.params['Wya'], self.params['by']
        Tx = X.shape[2]
        
        grads = {}
        grads['dWya'], grads['dby'] = np.zeros_like(Wya), np.zeros_like(by)
        grads['dWaa'], grads['da_prev'] = np.zeros_like(Waa), np.zeros_like(self.a0)
        grads['dWax'], grads['dba'] = np.zeros_like(Wax), np.zeros_like(ba)
        
        for t in reversed(range(Tx)):
            # compute gradient of loss wrt score
            dz = yhat[:, :, t] - y[:, :, t]
            cell_cache = X[:, :, t], a[:, :, t], a[:, :, t-1]
            grads = self.RNN_cell_backward(dz, grads, cell_cache)
        
        return grads, a
    
    
    def update_parameters(self, grads, lr):
        self.params['Wax'] -= lr * grads['dWax']
        self.params['Waa'] -= lr * grads['dWaa']
        self.params['Wya'] -= lr * grads['dWya']
        self.params['ba']  -= lr * grads['dba']
        self.params['by']  -= lr * grads['dby']
    
    
    def clip(self, grads, maxValue):
        for key in ['dWax', 'dWaa', 'dWya', 'dba', 'dby']:
            gradient = grads[key]
            grads[key] = np.clip(gradient, -maxValue, maxValue, out=gradient)
        return grads
    
    
    def train(self, X, y, lr, iters=1):
        loss_list = []
        for it in range(iters):
            loss, cache = self.RNN_forward(X, y)
            grads, a = self.RNN_backward(X, y, cache)
            # Clip gradients between -5 (min) and 5 (max)
            grads = self.clip(grads, 5)
            self.update_parameters(grads, lr)
            loss_list.append(loss)
        return loss, grads, a