Andrew Ng课程作业第三周(python版)

作业目的：学习针对multi-class classification的逻辑回归

 

作业内容：根据图片的400个像素点判定数字是几？（1-10）

20*20像素的灰度图片，每个值是像素点的亮度，共400个维度，值在0-255之间；如果用彩色图片，RGB，维度就是20*20*3=1200。

这类分类问题，boundary decision曲线必然由多项式组成，考虑所有的二次项，最终的多项式包含的二项式features有O(N*N)=N*N/2~400*400/2个，二项式的组成元素相当多。

提供的数据：ex3data1.txt（mooc可下载），数据集中有5000个samples，X是20*20个像素点组成的包含400个元素的矢量。y表征分类，数字从1-10,10用0表示

对假设函数/激活函数的理解：

方法：

处理过程：

方法一：

• step1:读取数据集

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.io import loadmat
%matplotlib inline

data = loadmat('d:\jupyter\ipython-notebooks-master\data\ex3data1.mat')
print(data)
data['X'].shape, data['y'].shape

((5000, 400), (5000, 1))

• step2:定义带正则化因子的cost function

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

def cost(theta, X, y, learningRate):
    theta = np.matrix(theta)
    X = np.matrix(X)
    y = np.matrix(y)
    first = np.multiply(-y, np.log(sigmoid(X * theta.T)))
    second = np.multiply((1 - y), np.log(1 - sigmoid(X * theta.T)))
    reg = (learningRate / 2 * len(X)) * np.sum(np.power(theta[:,1:theta.shape[1]], 2))
    return np.sum(first - second) / (len(X)) + reg

• step3:定义梯度函数

def gradient(theta, X, y, learningRate):
    theta = np.matrix(theta)
    X = np.matrix(X)
    y = np.matrix(y)
    
    parameters = int(theta.ravel().shape[1])
    error = sigmoid(X * theta.T) - y
    
    grad = ((X.T * error) / len(X)).T + ((learningRate / len(X)) * theta)
    
    # intercept gradient is not regularized
    grad[0, 0] = np.sum(np.multiply(error, X[:,0])) / len(X)
    
    return np.array(grad).ravel()

• step 4:one_vs_all，将y的1-10分类为 class i(1)和not class i(0)，选择是i的值为1，非i的值为0，这样得到h(x)的值越接近1的说明对应的场景下i的取值为输出的结果【max(h(x))】

from scipy.optimize import minimize
def one_vs_all(X, y, num_labels, learning_rate):
     rows = X.shape[0]      #sample数
     params= X.shape[1]   #X的参数个数/feature数
   
     #初始化theta，考虑theta(0)，所以theta的列数为params+1
     all_theta = np.zeros((num_labels, params+1))
  
     #X需要插入一列
     X = np.insert(X,0,values=np.ones(rows),axis=1)

     #y有1-10 10种不同的类别，one_vs_all分类将得到10种分类结果，theta有10组不同的结果
     # labels are 1-indexed instead of 0-indexed
     for i in range(1, num_labels + 1):
        theta = np.zeros(params + 1)
        y_i = np.array([1 if label == i else 0 for label in y])
        y_i = np.reshape(y_i, (rows, 1))
        # minimize the objective function
        fmin = minimize(fun=cost, x0=theta, args=(X, y_i, learning_rate), method='TNC', jac=gradient_with_loop)
        all_theta[i-1,:] = fmin.x
    
     return all_theta

• step5:初始化参数

rows = data['X'].shape[0]
params = data['X'].shape[1]

all_theta = np.zeros((10, params + 1))

X = np.insert(data['X'], 0, values=np.ones(rows), axis=1)

theta = np.zeros(params + 1)

#初始label=0的赋值1，其他为0
y_0 = np.array([1 if label == 0 else 0 for label in data['y']])
y_0 = np.reshape(y_0, (rows, 1))

X.shape, y_0.shape, theta.shape, all_theta.shape

• step6:得到theta值

all_theta = one_vs_all(data['X'], data['y'], 10, 1)
all_theta

• step7:预测

def predict_all(X, all_theta):
    rows = X.shape[0]
    params = X.shape[1]
    num_labels = all_theta.shape[0]
    
    # same as before, insert ones to match the shape
    X = np.insert(X, 0, values=np.ones(rows), axis=1)
    
    # convert to matrices
    X = np.matrix(X)
    all_theta = np.matrix(all_theta)
    
    # compute the class probability for each class on each training instance
    h = sigmoid(X * all_theta.T)
    print(h.shape)
    
    # create array of the index with the maximum probability,取每一行的最大值的索引
    h_argmax = np.argmax(h, axis=1)

    # because our array was zero-indexed we need to add one for the true label prediction，0-9改为1-10
    h_argmax = h_argmax + 1
    
    return h_argmax

• step8:计算精确度

y_pred = predict_all(data['X'], all_theta)
correct = [1 if a == b else 0 for (a, b) in zip(y_pred, data['y'])]
accuracy = (sum(map(int, correct)) / float(len(correct)))
print('accuracy = {0}%'.format(accuracy * 100))

方法二：

使用skicit-learn

from sklearn.model_selection import train_test_split
X_train,X_test, y_train, y_test = train_test_split(X, y, random_state=1) 
print(X_train.shape) 
print(y_train.shape)  
print(X_test.shape) 
print(y_test.shape)  
LogReg = LogisticRegression()
model=LogReg.fit(X_train,y_train)
y_pred=model.predict(X_test)
print(u"预测准确度为：%f%%"%np.mean(np.float64(y_pred == y_test)*100))