minist手写数据集的识别

#导入包
from mnist import MNIST
import matplotlib.pyplot as plt
import numpy as np
import random
from scipy.spatial.distance import cdist
import copy
from prettytable import PrettyTable
np.set_printoptions(precision=32)
#读取数据
mndata = MNIST('DataSet')
images, labels = mndata.load_training()
images_list=[[],[],[],[],[],[],[],[],[],[]]
print(len(images_list))
for image_id in range(len(images)):    
    images_list[labels[image_id]].append(images[image_id])
print(len(images_list[0]))
print(len(images_list[1]))
print(len(images_list[2]))
print(len(images_list[3]))
print(len(images_list[4]))
print(len(images_list[5]))
print(len(images_list[6]))
print(len(images_list[7]))
print(len(images_list[8]))
print(len(images_list[9]))
#以上根据label将sample分放在10个子列表中,每个子列表中存放对应数字类的所有样本。
##以下为your code，生成images_train，images_test两个四维数组，分别存放数字0-9的训练集以及测试集。images_train的尺寸是10,800,28,28.
##images_test的尺寸是10,200,28,28.
##1）推荐使用numpy函数random.choice或者random.sample来随机选取某类数字任意1000张中的800张组成训练集，剩余200张组成测试集。
##2）像素数值归一化到0-1之间
##此步骤 4分
images_choose_1000 = [[], [], [], [], [], [], [], [], [], []]
#从0-9十个数字中选取1000张图像存放
for i in range(10):
    images_choose_1000[i] = random.sample(images_list[i], 1000)
#选取800个作为训练集，200个作为测试集
index_800 = [[], [], [], [], [], [], [], [], [], []]
index_200 = [[], [], [], [], [], [], [], [], [], []]
images_train = [[], [], [], [], [], [], [], [], [], []]
images_test = [[], [], [], [], [], [], [], [], [], []]
a = np.arange(1000)
for i in range(10):
    index_800[i] = np.random.choice(a, size = 800, replace = False) #存取训练集800个下标
    index_200[i] = list(set(a) - set(index_800[i]))  #获取其余200个下标
    for j in index_800[i]:
        images_train[i].append(images_choose_1000[i][j])
    for j in index_200[i]:
        images_test[i].append(images_choose_1000[i][j])
images_train = np.reshape(np.asarray(images_train), (10,-1,28,28)) / 255 #归一化，iamges_train 存放训练数据
images_test = np.reshape(np.asarray(images_test), (10,-1,28,28)) / 255  #images_test存放测试数据
# 计算并画出数字“3”training set以及testing set的均值图像。 
#此步骤得2分
print(images_train.shape)
print(images_test.shape)
nub=3 #数字3
##以下your code。可使用plt.imshow画图。
images_train_mean = np.zeros((10,28, 28)) #images_train_mean --- 10个数字的均值图像矩阵
images_test_mean = np.zeros((10,28, 28))  
images_train_mean = np.mean(images_train, axis=1)
images_test_mean = np.mean(images_test,axis= 1)
plt.figure()
plt.subplot(1,2,1)
plt.imshow(images_train_mean[nub], cmap = 'gray')
plt.title("Mean_value Fig-traing_set_3")
plt.subplot(1,2,2)
plt.imshow(images_test_mean[nub], cmap = 'gray')
plt.title("Mean_value Fig-test_set_3")
plt.show()
positive_nub=3 #以positive_nub为images_train第一轴的索引，获得“3”为正样本，以negative_nub为索引获得数字“6”为负样本
negative_nub=6
negative_train_sample = np.ones((800, 785, 1))
##以下your code，获得投影方向W之后，将共1600个正、负样本投影到该方向，并统计正样本投影结果的分布以及负样本投影结果的分布，绘制在一张图中，由此寻找恰当的threshold。
##最后，调用images_test中正负样本集，通过W'X是否大于threshold来判断正负性，并统计“3”“6”的分类error。请输出 confusion matrix。至此，任务一结束。
#此步骤得4分。

positive_mean = np.zeros((785, 1))
negative_mean = np.zeros((785, 1))
positive_mean = np.mean(positive_train_sample, axis=0)  #positive_mean存取数字3的均值向量
negative_mean = np.mean(negative_train_sample, axis=0)  #negative_mean存取数字6的均值向量

#求SB --- 类间散度矩阵
SB = np.dot(np.subtract(negative_mean,positive_mean), np.subtract(negative_mean, positive_mean).T)
#求SW --- 类内散度矩阵
SW = np.zeros((785, 785))
for i in range(800):
   SW = np.add(SW + np.dot((positive_train_sample[i] - positive_mean), (positive_train_sample[i] - positive_mean).T ) , np.dot((negative_train_sample[i] - negative_mean), (negative_train_sample[i] - negative_mean).T ))
#使用奇异值分解求解SW的逆矩阵
U, sigma, VT = np.linalg.svd(SW)
sigma_diag = np.diag(sigma)
SW_inv = np.dot(np.dot((VT).T, np.linalg.inv(sigma_diag)), (U).T) #SW_inv即SW的逆矩阵
#求最佳的投影方向w
w = np.dot(SW_inv, np.subtract(negative_mean, positive_mean)) 
positive_result = np.zeros((800, 1))# 3投影
negative_result = np.zeros((800, 1)) # 6 投影
for i in range(800):
    positive_result[i][0] = np.dot(w.T, positive_train_sample[i])
    negative_result[i][0] = np.dot(w.T, negative_train_sample[i])
plt.hist(positive_result, bins = 13, facecolor="blue", edgecolor="black", label = "Number 3", density = "True") #蓝色的是数字3
plt.hist(negative_result, bins = 13, facecolor="red", edgecolor="black",  label = "Number 6", density = "True") #红色的是数字6
plt.text(-0.02,80,"Number 3")
plt.text(0.02,75, "Number 6")
# 显示横轴标签
plt.xlabel("The range of x")
# 显示纵轴标签
plt.ylabel("The frequency of x")
# 显示图标题
plt.title("Projection distribution histogram of numbers 3 and 6")
plt.show()
#如何确定阈值---选取数字3投影后取值最大的并总数的5%和数字6投影后取值最小的并占总数的5%
#取这些数字的平均值作为Threshold
positive_result_sort = np.sort(positive_result[:, 0])
negative_result_sort = np.sort(negative_result[:, 0])  #对正类和负类进行排序
positive_result_sort_inverse = positive_result_sort[::-1] #对正类进行倒置
Threshold = 0.0
for i in range(int(800*0.05)):
    Threshold = Threshold + positive_result_sort_inverse[i] + negative_result_sort[i]
Threshold = Threshold / (2 * 800 * 0.05)  #获得阈值

positive_test_sample_temp = images_test[3].reshape(200, 784, 1)
negative_test_sample_temp = images_test[6].reshape(200, 784, 1)

##每个xi最后面应该加1
positive_test_sample = np.ones((200, 785, 1))
negative_test_sample = np.ones((200, 785, 1))
for i in range(200):
    for j in range(784):
        positive_test_sample[i][j][0] = positive_test_sample_temp[i][j][0]
        negative_test_sample[i][j][0] = negative_test_sample_temp[i][j][0]
##开始训练
TP = 0
FN = 0
FP = 0
TN = 0
for i in range(200):
    if np.dot(w.T, positive_test_sample[i]) < Threshold :
        TP = TP + 1
    else:
        FN = FN + 1
    if np.dot(w.T, negative_test_sample[i]) > Threshold:
        TN = TN + 1
    else:
        FP = FP + 1
Error_rate = (FP + FN ) / 400
##以下任务二的your code。请输出三维空间的数字投影图，以不同颜色区分数字。 请在W1，W2， W3确定的三维空间中，计算“3”，“6”，以及“未知数字”三类之间Sw，Sb，
##以数字类为横轴，以Sb/Sw为纵轴，绘制曲线，并确定能与“3”“6”最大程度区分的数字。 
##此步骤得6分
images_choose_5000 = [[], [], [], [], [], [], [], [], [], []]  #images_choose_5000 每个数字取5000张
for i in range(10):
    images_choose_5000[i] = random.sample(images_list[i], 5000)
images_choose_5000_array = np.reshape(images_choose_5000, (10, -1, 784, 1)) / 255 #归一化
images_choose_5000_sample = np.ones((10, 5000, 785, 1))
for i in range(10):
    for j in range(5000):
        for k in range(784):
            images_choose_5000_sample[i][j][k][0] = images_choose_5000_array[i][j][k][0]  #对所获得的数据升维操作（加1操作）
#计算各类数字样本的类内散度Swi
#1,先计算各类数字的均值矩阵 并画出各个均值图像
mean_array = np.zeros((10,785,1))  
for i in range(10):
    mean_array[i] = np.mean(images_choose_5000_sample[i], axis = 0) #mean_array存取每个数字的均值矩阵
plt.figure()
for i in range(10):
    plt.subplot(2,5,i+1)
    plt.imshow(np.reshape(mean_array[i][0:784], (28, 28)), cmap = 'gray') #画出各个数字的均值图像
#2,计算各类数字的内散度矩阵 --- SW_array存储的是每个数字的内散度矩阵
SW_array = np.zeros((10, 785, 785))
for i in range(10):
    for j in range(5000):
        SW_array[i] = SW_array[i] + np.dot((images_choose_5000_sample[i][j] - mean_array[i]), (images_choose_5000_sample[i][j] - mean_array[i]).T)
#SW_sum为每个数字对应的类内散度矩阵之和 --- 全局类内散度
SW_sum = np.zeros((785,785))  
for i in range(10):
    SW_sum = SW_sum + SW_array[i]
    #计算Sb_task2 --- 任务的2的类间矩阵
Sb_task2 = np.zeros((785, 785)) #初始化
mean_all = np.zeros((785, 1)) 
mean_all = np.mean(mean_array, axis = 0) #10数字合并的均值矩阵
for i in range(10):
    Sb_task2 = Sb_task2 + 5000 * np.dot((mean_array[i] - mean_all), (mean_array[i] - mean_all).T) #算出十分类的类间散度矩阵
    #下面求SW_sum的逆矩阵,使用奇异值分解
U_task2, sigma_task2, VT_task2 = np.linalg.svd(SW_sum)   
sigma_task2_diag = np.diag(sigma_task2)
SW_sum_inv = np.matmul(np.matmul((VT_task2).T, np.linalg.inv(sigma_task2_diag)), (U_task2).T) #SW_sum_inv 全局类间散度矩阵的逆矩阵
#下面求SW_sum_inv * Sb_task2的特征值和特征向量
Feature_value = np.zeros((785, 1)) #Feature_value --- 特征值
Feature_vector = np.zeros((785, 1))#Feature_vector --- 特征向量
Feature_value, Feature_vector = np.linalg.eig(np.matmul(SW_sum_inv, Sb_task2))
#在pycharm观察到最大的三个特征值（实数）对应的特征向量也是实数部分
Feature_value_real = np.real(Feature_value) #取特征值的实数部分
Feature_vector_real = np.real(Feature_vector) #取特征限量的实数部分
index_list = np.argsort(Feature_value)
max_list_value = np.asarray([Feature_value_real[index_list[-1]], Feature_value_real[index_list[-2]], Feature_value_real[index_list[-3]]])
max_feature_vector = np.zeros((785,3))
max_feature_vector[:,0] = Feature_vector_real[:, index_list[-1]]
max_feature_vector[:,1] = Feature_vector_real[:, index_list[-2]]
max_feature_vector[:,2] = Feature_vector_real[:, index_list[-3]] #max_feature_vector存储的就是所取得三个特征向量（已验证）
#将10个数字全部投影到这个三维空间中
abb_array = np.zeros((10, 5000, 3)) #abb_array存储每个数字对应在三维空间的坐标
for i in range(10):
    for j in range(5000):
        for k in range(3):
            abb_array[i][j][k] = np.dot(max_feature_vector[:,k].T, images_choose_5000_sample[i][j]) #投影降维
from mpl_toolkits.mplot3d import Axes3D
ax = plt.subplot(projection = '3d')
for i in range(10): #画出10个数字在降维空间内的投影图
    ax.scatter(abb_array[i][:, 0], abb_array[i][:, 1], abb_array[i][:, 2], s = np.ones(5000)/4, label = str(i))
ax.set_xlabel('projection value in w1 direction')
ax.set_ylabel('projection value in w2 direction')
ax.set_zlabel('projection value in w3 direction')
ax.set_title("projection of 10 digits")
plt.legend(loc = 'upper left')
#计算8个数字和3，6数字之间的类间散度，mean_array (10,785,1)保存了各个数字的均值限量
Sb_array_3 = np.zeros((10, 785, 785))   #Sb_array 存储10个数字三分类问题的Sb矩阵（多了3,3,6和3，6,6）
mean_3_array = np.zeros((10, 785, 1))
#1 下面开始算每个三分类矩阵的均值矩阵
for i in range(10):
    mean_3_array[i] = (mean_array[3] + mean_array[6] + mean_array[i]) / 3 #mean_3_array每个三分类的均值矩阵
for i in range(10):
        Sb_array_3[i] = Sb_array_3[i] + 5000 * np.dot((mean_array[3] - mean_3_array[i]),(mean_array[3] - mean_3_array[i]).T) + 5000*np.dot((mean_array[6] - mean_3_array[i]),(mean_array[6] - mean_3_array[i]).T) + 5000 * np.dot((mean_array[i] - mean_3_array[i]),(mean_array[i] - mean_3_array[i]).T)
#2 下面开始计算Sw_array_3
Sw_array_3 = np.zeros((10,785, 785)) #Sw_array_3存储的十个三分类的Sw
for i in range(10):
    Sw_array_3[i] = SW_array[i] + SW_array[3] + SW_array[6]
index_list_3 = [0, 1, 2, 4, 5, 7, 8, 9]
tr_up = np.zeros((10, 1))
tr_bottom = np.zeros((10, 1))
result_extra_number = np.zeros((10, 1))
for i in range(10):
    for j in range(3):
        tr_up[i] = tr_up[i] + np.dot(np.dot(max_feature_vector[:, j].T, Sb_array_3[i]), max_feature_vector[:, j])
        tr_bottom[i] = tr_bottom[i] + np.dot(np.dot(max_feature_vector[:, j].T, Sw_array_3[i]), max_feature_vector[:, j])
    result_extra_number[i] = tr_up[i] / tr_bottom[i]
 #将result_extra_number中除去3和6的J值存入result_extra_number_final中
result_extra_number_final = []
for i in range(10):
    if i == 3 or i == 6:
        continue
    else:
        result_extra_number_final.append(result_extra_number[i][0])
#下面开始画曲线图
plt.figure()
plt.plot(index_list_3, result_extra_number_final)
new_ticks = np.linspace(0, 9, 10)
plt.xticks(new_ticks)
plt.xlabel("Type of number")
plt.ylabel("Sb/Sw")
plt.title("Sb/Sw value of each number except 3 and 6")
plt.show()