题目：
要求：在UCI数据集上的Iris和sonar数据上验证算法的有效性；Iris数据3类，4维，150个数据；Sonar数据2类，60维，208个样本；
训练和测试样本有三种方式进行划分：（三选一）
1）将数据随机分训练和测试，多次平均求结果
2）k折交叉验证
3）留1法
（针对不同维数，画出曲线图；）
目录
一、基础知识
二、仿真实验
基础知识
如图所示，Fisher线性判别分析的最重要的作用是将数据降维。用一句话总结就是：用数学上的闭合解，去二分类问题。所以Fisher的核心问题是，这个确切的闭合解是如何求得的。首先引入几个必要的概念：均值、类内离散度、类间离散度。
作为二分类问题，期望的是类内离散度最小，类间离散度最大。所以准则函数可以定义如下：
对于J(w)的求解可以概括为两步：1.将J(w)函数转换为显性方程。2.使用拉格朗日乘子法求梯度。
通过计算可以得知，w*只和Sw^-1*(u1 – u2)相关。
至此就可以开始写程序了。
仿真实验
2.1对数据集的分析
数据集采用sonar以及iris两个数据集。
Sonar：两类； 208行61列，对应前97行R（Rock）和 后111行M（Mine）。
Iris：三类； 150行5列，前4列分别为花萼长度、花萼宽度、花瓣长度、花瓣宽度； 最后一列为标签。每类50个样本。
本实验中两数据集上都采用留一法验证。
对于Sonar数据集来说，由于sonar一共有60维，因此程序中分别取1，2，3…60维，分别验证每个维度对正确率的影响。
对于iris数据集来说，由于Fisher只可以进行二分类问题，因此将数据分为了12，13，23三个二分类问题。
2.2 Sonar代码
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
def Fisher(x1, x2, n, c):
#计算各类均值
u1 = np.mean(x1, axis=0)
u2 = np.mean(x2, axis=0)
# u1 = u1.reshape(-1, 1)
# u2 = u2.reshape(-1, 1)
#计算S1，S2,类内离散度矩阵
S1 = np.zeros((n, n))
S2 = np.zeros((n, n))
if c == 0: # 第一种情况
for i in range(0,96):
S1 += (x1[i].reshape(-1, 1)-u1).dot((x1[i].reshape(-1, 1)-u1).T)
for i in range(0,111):
S2 += (x2[i].reshape(-1, 1)-u2).dot((x2[i].reshape(-1, 1)-u2).T)
if c == 1:
for i in range(0,97):
S1 += (x1[i].reshape(-1, 1)-u1).dot((x1[i].reshape(-1, 1)-u1).T)
for i in range(0,110):
S2 += (x2[i].reshape(-1, 1)-u2).dot((x2[i].reshape(-1, 1)-u2).T)
#计算Sw
Sw = S1 + S2
#计算W以及W0
W = np.linalg.inv(Sw) @ (u1 - u2)
u_1 = u1.T @ W #投影在一维的均值
u_2 = u2.T @ W
W0 = -0.5 * (u_1 + u_2) #分界点
return W, W0
def Classify(W, W0, test):
y = W0 + test @ W
if y > 0:
return 1
else:
return 0
sonar = pd.read_csv("你的文件地址",header=None,sep=',')
sonar1 = sonar.iloc[:, 0:60]
sonar2 = np.mat(sonar1)
ACC = np.zeros(60)
for n in range(1, 61):
sonar_random = (np.random.permutation(sonar2.T)).T
P1 = sonar_random[0:97, 0:n]
P2 = sonar_random[97:208, 0:n]
acc = np.zeros(10)
for t in range(10):
count = 0
for i in range(208):
if i < 97:
test = P1[i]
P1_real = np.delete(P1, i, axis=0)
W, W0 = Fisher(P1_real, P2, n, 0)
if Classify(W, W0, test):
count += 1
else:
test = P2[i-97]
P2_real = np.delete(P2, i-97, axis=0)
W, W0 = Fisher(P1, P2_real, n, 1)
if not Classify(W, W0, test):
count += 1
acc[t] = count / 208
ACC[n-1] = acc.mean()
print("当前为", n, "维，正确率是", ACC[n-1] )
x = np.arange(1, 61, 1)
plt.xlabel("Dim")
plt.ylabel("Accuracy")
plt.plot(x, ACC, 'r')
plt.savefig('sonar-result.png',dpi=2000)
2.2 iris代码
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
def Fisher(X1,X2,n,c):
# 计算三类样本的类均值向量
m1=(np.mean(X1,axis = 0))
m2=(np.mean(X2,axis = 0))
m1 = m1.reshape(n,1) # 将行向量转换为列向量以便于计算
m2 = m2.reshape(n,1)
#计算类内离散度矩阵
S1 = np.zeros((n,n)) # m1 = within_class_scatter_matrix1
S2 = np.zeros((n,n)) # m2 = within_class_scatter_matrix2
if c == 0: # 第一种情况
for i in range(0,49):
S1 += (X1[i].reshape(n,1)-m1).dot((X1[i].reshape(n,1)-m1).T)
for i in range(0,50):
S2 += (X2[i].reshape(n,1)-m2).dot((X2[i].reshape(n,1)-m2).T)
if c == 1:
for i in range(0,50):
S1 += (X1[i].reshape(n,1)-m1).dot((X1[i].reshape(n,1)-m1).T)
for i in range(0,49):
S2 += (X2[i].reshape(n,1)-m2).dot((X2[i].reshape(n,1)-m2).T)
#计算总类内离散度矩阵S_w
S_w = S1 + S2
#计算最优投影方向 W
W = np.linalg.inv(S_w).dot(m1 - m2)
#在投影后的一维空间求两类的均值
m_1 = (W.T).dot(m1)
m_2 = (W.T).dot(m2)
#计算分类阈值 W0(为一个列向量)
W0 = -0.5*(m_1 + m_2)
return W,W0
def Classify(X,W,W0):
y = (W.T).dot(X) + W0
return y
#导入sonar.all-data数据集
iris = pd.read_csv('你的代码',header=None,sep=',')
iris1 = iris.iloc[0:150,0:4]
iris2 = np.mat(iris1)
Accuracy = 0
accuracy_ = np.zeros(10)
P1 = iris2[0:50,0:4]
P2 = iris2[50:100,0:4]
P3 = iris2[100:150,0:4]
G121 = np.ones(50)
G122 = np.ones(50)
G131 = np.zeros(50)
G132 = np.zeros(50)
G231 = np.zeros(50)
G232 = np.zeros(50)
# 留一法验证准确性
# 第一类和第二类的线性判别
count = 0
for i in range(100):
if i <= 49:
test = P1[i]
test = test.reshape(4,1)
train = np.delete(P1,i,axis=0) # 训练样本是一个列数为t的矩阵
W,W0 = Fisher(train,P2,4,0)
if (Classify(test,W,W0)) >= 0:
count += 1
G121[i] = Classify(test,W,W0)
else:
test = P2[i-50]
test = test.reshape(4,1)
train = np.delete(P2,i-50,axis=0)
W,W0 = Fisher(P1,train,4,1)
if (Classify(test,W,W0)) < 0:
count += 1
G122[i-50] = Classify(test,W,W0)
Accuracy12 = count/100
print("第一类和二类的分类准确率为:%.3f"%(Accuracy12))
# 第一类和第三类的线性判别
count = 0
for i in range(100):
if i <= 49:
test = P1[i]
test = test.reshape(4,1)
train = np.delete(P1,i,axis=0) # 训练样本是一个列数为t的矩阵
W,W0 = Fisher(train,P3,4,0)
if (Classify(test,W,W0)) >= 0:
count += 1
G131[i] = Classify(test,W,W0)
else:
test = P3[i-50]
test = test.reshape(4,1)
train = np.delete(P3,i-50,axis=0)
W,W0 = Fisher(P1,train,4,1)
if (Classify(test,W,W0)) < 0:
count += 1
G132[i-50] = Classify(test,W,W0)
Accuracy13 = count/100
print("第一类和三类的分类准确率为:%.3f"%(Accuracy13))
# 第二类和第三类的线性判别
count = 0
for i in range(100):
if i <= 49:
test = P2[i]
test = test.reshape(4,1)
train = np.delete(P2,i,axis=0) # 训练样本是一个列数为t的矩阵
W,W0 = Fisher(train,P3,4,0)
if (Classify(test,W,W0)) >= 0:
count += 1
G231[i] = Classify(test,W,W0)
else:
test = P3[i-50]
test = test.reshape(4,1)
train = np.delete(P3,i-50,axis=0)
W,W0 = Fisher(P2,train,4,1)
if (Classify(test,W,W0)) < 0:
count += 1
G232[i-50] = Classify(test,W,W0)
Accuracy23 = count/100
print("第二类和三类的分类准确率为:%.3f"%(Accuracy23))
# 画相关的图
import matplotlib.pyplot as plt
y1 = np.zeros(50)
y2 = np.zeros(50)
plt.figure(1)
plt.ylim((-0.5,0.5)) # y坐标的范围
#画散点图
plt.scatter(G121, y1,c='red', alpha=1, marker='.')
plt.scatter(G122, y2,c='k', alpha=1, marker='.')
plt.xlabel('Class:1-2')
plt.savefig('iris 1-2.png',dpi=2000)
plt.figure(2)
plt.ylim((-0.5,0.5)) # y坐标的范围
#画散点图
plt.scatter(G131, y1,c='red', alpha=1, marker='.')
plt.scatter(G132, y2,c='k', alpha=1, marker='.')
plt.xlabel('Class:1-3')
plt.savefig('iris 1-3.png',dpi=2000)
plt.figure(3)
plt.ylim((-0.5,0.5)) # y坐标的范围
#画散点图
plt.scatter(G231, y1,c='red', alpha=1, marker='.')
plt.scatter(G232, y2,c='k', alpha=1, marker='.')
plt.xlabel('Class:2-3')
plt.savefig('iris 2-3.png',dpi=2000)
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
def Fisher(X1,X2,n,c):
# 计算三类样本的类均值向量
m1=(np.mean(X1,axis = 0))
m2=(np.mean(X2,axis = 0))
m1 = m1.reshape(n,1) # 将行向量转换为列向量以便于计算
m2 = m2.reshape(n,1)
#计算类内离散度矩阵
S1 = np.zeros((n,n)) # m1 = within_class_scatter_matrix1
S2 = np.zeros((n,n)) # m2 = within_class_scatter_matrix2
if c == 0: # 第一种情况
for i in range(0,49):
S1 += (X1[i].reshape(n,1)-m1).dot((X1[i].reshape(n,1)-m1).T)
for i in range(0,50):
S2 += (X2[i].reshape(n,1)-m2).dot((X2[i].reshape(n,1)-m2).T)
if c == 1:
for i in range(0,50):
S1 += (X1[i].reshape(n,1)-m1).dot((X1[i].reshape(n,1)-m1).T)
for i in range(0,49):
S2 += (X2[i].reshape(n,1)-m2).dot((X2[i].reshape(n,1)-m2).T)
#计算总类内离散度矩阵S_w
S_w = S1 + S2
#计算最优投影方向 W
W = np.linalg.inv(S_w).dot(m1 - m2)
#在投影后的一维空间求两类的均值
m_1 = (W.T).dot(m1)
m_2 = (W.T).dot(m2)
#计算分类阈值 W0(为一个列向量)
W0 = -0.5*(m_1 + m_2)
return W,W0
def Classify(X,W,W0):
y = (W.T).dot(X) + W0
return y
#导入sonar.all-data数据集
iris = pd.read_csv('F:\\1_数模\\PY_CODE\\Fisher线性判别分析\\iris\\iris.data',header=None,sep=',')
iris1 = iris.iloc[0:150,0:4]
iris2 = np.mat(iris1)
Accuracy = 0
accuracy_ = np.zeros(10)
P1 = iris2[0:50,0:4]
P2 = iris2[50:100,0:4]
P3 = iris2[100:150,0:4]
G121 = np.ones(50)
G122 = np.ones(50)
G131 = np.zeros(50)
G132 = np.zeros(50)
G231 = np.zeros(50)
G232 = np.zeros(50)
# 留一法验证准确性
# 第一类和第二类的线性判别
count = 0
for i in range(100):
if i <= 49:
test = P1[i]
test = test.reshape(4,1)
train = np.delete(P1,i,axis=0) # 训练样本是一个列数为t的矩阵
W,W0 = Fisher(train,P2,4,0)
if (Classify(test,W,W0)) >= 0:
count += 1
G121[i] = Classify(test,W,W0)
else:
test = P2[i-50]
test = test.reshape(4,1)
train = np.delete(P2,i-50,axis=0)
W,W0 = Fisher(P1,train,4,1)
if (Classify(test,W,W0)) < 0:
count += 1
G122[i-50] = Classify(test,W,W0)
Accuracy12 = count/100
print("第一类和二类的分类准确率为:%.3f"%(Accuracy12))
# 第一类和第三类的线性判别
count = 0
for i in range(100):
if i <= 49:
test = P1[i]
test = test.reshape(4,1)
train = np.delete(P1,i,axis=0) # 训练样本是一个列数为t的矩阵
W,W0 = Fisher(train,P3,4,0)
if (Classify(test,W,W0)) >= 0:
count += 1
G131[i] = Classify(test,W,W0)
else:
test = P3[i-50]
test = test.reshape(4,1)
train = np.delete(P3,i-50,axis=0)
W,W0 = Fisher(P1,train,4,1)
if (Classify(test,W,W0)) < 0:
count += 1
G132[i-50] = Classify(test,W,W0)
Accuracy13 = count/100
print("第一类和三类的分类准确率为:%.3f"%(Accuracy13))
# 第二类和第三类的线性判别
count = 0
for i in range(100):
if i <= 49:
test = P2[i]
test = test.reshape(4,1)
train = np.delete(P2,i,axis=0) # 训练样本是一个列数为t的矩阵
W,W0 = Fisher(train,P3,4,0)
if (Classify(test,W,W0)) >= 0:
count += 1
G231[i] = Classify(test,W,W0)
else:
test = P3[i-50]
test = test.reshape(4,1)
train = np.delete(P3,i-50,axis=0)
W,W0 = Fisher(P2,train,4,1)
if (Classify(test,W,W0)) < 0:
count += 1
G232[i-50] = Classify(test,W,W0)
Accuracy23 = count/100
print("第二类和三类的分类准确率为:%.3f"%(Accuracy23))
# 画相关的图
import matplotlib.pyplot as plt
y1 = np.zeros(50)
y2 = np.zeros(50)
plt.figure(1)
plt.ylim((-0.5,0.5)) # y坐标的范围
#画散点图
plt.scatter(G121, y1,c='red', alpha=1, marker='.')
plt.scatter(G122, y2,c='k', alpha=1, marker='.')
plt.xlabel('Class:1-2')
plt.savefig('iris 1-2.png',dpi=2000)
plt.figure(2)
plt.ylim((-0.5,0.5)) # y坐标的范围
#画散点图
plt.scatter(G131, y1,c='red', alpha=1, marker='.')
plt.scatter(G132, y2,c='k', alpha=1, marker='.')
plt.xlabel('Class:1-3')
plt.savefig('iris 1-3.png',dpi=2000)
plt.figure(3)
plt.ylim((-0.5,0.5)) # y坐标的范围
#画散点图
plt.scatter(G231, y1,c='red', alpha=1, marker='.')
plt.scatter(G232, y2,c='k', alpha=1, marker='.')
plt.xlabel('Class:2-3')
plt.savefig('iris 2-3.png',dpi=2000)
三、实验结果分析
3.1Sonar
实验中，每个维度上都对数据集的正确率进行了10次求平均，保证数据的稳定性。并最终把每个维度的平均正确率存在了数组ACC中。最终以ACC为x轴，准确率为y轴，进行绘图。
由图可见，可以大致的认为准确率和维度成正比关系。
3.2 iris
在iris数据集的实验中，主要对12，13，23类别进行了划分。并将每次划分后的结果在一维上，以红色和黑色的点呈现出来。
实验结果如上图所示，12，13类可以非常完美的分开；而23类分类时，两类降维后依然挨得很近，所以准确率略低。
四、总结及疑问
在此次实验中，感受到了使用闭合解的优势，求解过程的代码量极其的少。
但还是有些疑问，在求解Sw时，为什么必须将均值u1,u2以及样本x1,x2转换为列向量。在不转换为列向量时，使用np.linalg.inv()函数求解Sw的逆时会出现报错。
6025
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
