首先,绘制原始数据:
其中
X1=load('ex3red.dat');
X2=load('ex3green.dat');
X3=load('ex3blue.dat');
hold on
scatter(X1(:,1),X1(:,2),'r')
scatter(X2(:,1),X2(:,2),'g')
scatter(X3(:,1),X3(:,2),'blue')
grid on
mu1=mean(X1)
mu2=mean(X2)
mu3=mean(X3)
Sw=(X1-mu1)'*(X1-mu1)+(X3-mu3)'*(X3-mu3)
SW=(X1-mu1)'*(X1-mu1)+(X2-mu2)'*(X3-mu3)...
+(X3-mu3)'*(X3-mu3)
Sb=(mu1-mu3)'*(mu1-mu3)
SB=(mu1-mu3)'*(mu1-mu3)+(mu1-mu2)'*(mu1-mu2)...
+(mu2-mu3)'*(mu2-mu3)
其中
Sw = 2×2
4.8008 -1.4172
-1.4172 29.0730
SW = 2×2
6.2625 -2.8762
-1.7051 29.6357
Sb = 2×2
1.5482 3.2831
3.2831 6.9621
SB = 2×2
16.8330 18.0347
18.0347 22.2876
计算theta:theta=(Sw)\(mu1-mu3)'
此处没有使用inv(Sw),因为这样相比于直接反除会损失精度。
theta = 2×1
0.2901
0.1049
(3) 绘制图像如下:
其中,直线的斜率为theta(2)/theta(1)
根据两向量之间的夹角以及斜率计算投影之后的x,y坐标
计算方式如下:
%投影在直线上的长度
L1=(X1*theta)/sqrt(theta(2)^2+theta(1)^2)
L2=(X3*theta)/sqrt(theta(2)^2+theta(1)^2)
%斜率
k=theta(2)/theta(1)
x1=L1/sqrt(1+k^2)
x2=L2/sqrt(1+k^2)
y1=x1*k
y2=x2*k
绘制图像如下:
步骤如下:
此时计算了矩阵的特征值和特征向量,调用MATLAB中的库函数
[V,D]=eig(SW\SB) %其中 V是特征向量矩阵,D是特征值对角矩阵
w1=V(:,1)
w2=V(:,2)
其中w1和w2是两个特征向量。
V = 2×2
0.9670 -0.7437
0.2546 0.6685
D = 2×2
3.9201 0
0 0.0705
w1 = 2×1
0.9670
0.2546
w2 = 2×1
-0.7437
0.6685
对于第一个特征向量绘制图像如下:
对于第二个特征向量w2:
可以看到,显然我们要使用对应特征值最大的那个向量,即w1。
X1=load('ex3red.dat');
X2=load('ex3green.dat');
X3=load('ex3blue.dat');
hold on
scatter(X1(:,1),X1(:,2),'r')
scatter(X2(:,1),X2(:,2),'g')
scatter(X3(:,1),X3(:,2),'blue')
grid on
%two classes
figure
hold on
scatter(X1(:,1),X1(:,2),'r','filled')
scatter(X3(:,1),X3(:,2),'blue','filled')
title('red and blue')
xlabel('x')
ylabel('y')
xlim([0.00 10.00])
ylim([0.00 10.00])
mu1=mean(X1)
mu2=mean(X2)
mu3=mean(X3)
Sw=(X1-mu1)'*(X1-mu1)+(X3-mu3)'*(X3-mu3)
SW=(X1-mu1)'*(X1-mu1)+(X2-mu2)'*(X3-mu3)...
+(X3-mu3)'*(X3-mu3)
Sb=(mu1-mu3)'*(mu1-mu3)
SB=(mu1-mu3)'*(mu1-mu3)+(mu1-mu2)'*(mu1-mu2)...
+(mu2-mu3)'*(mu2-mu3)
[V,D]=eig(SW\SB) %其中 V是特征向量矩阵,D是特征值对角矩阵
w1=V(:,1)
w2=V(:,2)
theta=(Sw)\(mu1-mu3)'
disp(theta(1,1))
figure
hold on
scatter(X1(:,1),X1(:,2),'r','filled')
scatter(X3(:,1),X3(:,2),'blue','filled')
x=linspace(0,10,100);
y=(theta(2)/theta(1))*x
plot(x,y,'black')
title('LDA for two-classes')
xlabel('x')
ylabel('y')
xlim([0.00 10.00])
2-ylim([0.00 10.00])
%投影在直线上的长度
L1=(X1*theta)/sqrt(theta(2)^2+theta(1)^2)
L2=(X3*theta)/sqrt(theta(2)^2+theta(1)^2)
%斜率
k=theta(2)/theta(1)
x1=L1/sqrt(1+k^2)
x2=L2/sqrt(1+k^2)
y1=x1*k
y2=x2*k
figure
hold on
scatter(X1(:,1),X1(:,2),'r','filled')
scatter(X3(:,1),X3(:,2),'blue','filled')
scatter(x1',y1','red','x')
scatter(x2',y2','blue','x')
x=linspace(0,10,100);
y=(theta(2)/theta(1))*x
plot(x,y,'black')
title('LDA')
xlabel('x')
ylabel('y')
xlim([0.00 10.00])
ylim([0.00 10.00])
D1=(X1*w1)/sqrt(w1(2)^2+w1(1)^2)
D2=(X2*w1)/sqrt(w1(2)^2+w1(1)^2)
D3=(X3*w1)/sqrt(w1(2)^2+w1(1)^2)
K=w1(2)/w1(1)
x1=D1/sqrt(1+K^2)
x2=D2/sqrt(1+K^2)
x3=D3/sqrt(1+K^2)
y1=x1*K
y2=x2*K
y3=x3*K
figure
hold on
scatter(X1(:,1),X1(:,2),'r','filled')
scatter(X2(:,1),X2(:,2),'g','filled')
scatter(X3(:,1),X3(:,2),'blue','filled')
scatter(x1',y1','red','x')
scatter(x2',y2','green','x')
scatter(x3',y3','blue','x')
x=linspace(0,10,100);
y=(w1(2)/w1(1))*x
plot(x,y,'black')
D1=(X1*w2)/sqrt(w2(2)^2+w2(1)^2)
D2=(X2*w2)/sqrt(w2(2)^2+w2(1)^2)
D3=(X3*w2)/sqrt(w2(2)^2+w2(1)^2)
K=w2(2)/w2(1)
x1=D1/sqrt(1+K^2)
x2=D2/sqrt(1+K^2)
x3=D3/sqrt(1+K^2)
y1=x1*K
y2=x2*K
y3=x3*K
figure
hold on
scatter(X1(:,1),X1(:,2),'r','filled')
scatter(X2(:,1),X2(:,2),'g','filled')
scatter(X3(:,1),X3(:,2),'blue','filled')
scatter(x1',y1','red','x')
scatter(x2',y2','green','x')
scatter(x3',y3','blue','x')
x=linspace(0,10,100);
y=(w2(2)/w2(1))*x
plot(x,y,'black')