在 matlab 中绘制 3D 二次判别分析的输出

Question

我正在尝试在 MATLAB 中绘制二次判别分析的输出。我之前问过线性分析（LDA问题供参考），但是，现在我想使用二次判别分析进行研究。对于 LDA，我使用 fimplicit3 来绘制平原，但这需要知道要评估的函数的方程形式。我正在努力弄清楚这应该是什么。下面的代码适用于 LDA，但是当我更改为 QDA 时，它显然是不正确的（运行代码时您会看到，平原与数据相交），因为它不包括 QDA 引入的额外“二次系数”。我尝试了一些方法，但都没有成功，为了清楚起见，我省略了这些方法。

总而言之，我需要传递给第 135、137、139 行的“fimplicit3”以绘制弯曲平面的方程是什么？有更好的方法吗？

在代码示例下面，我包含了二维线性和二次方程，以防有帮助。

clear
clc
close all

%----------------------------
x1 = 1:0.01:3;
r = -1 + (1+1)*rand(1,201);
x1 = x1 + r;

y1 = 1:0.01:3;
r = -1 + (1+1)*rand(1,201);
y1 = y1 + r;

z1 = 1:0.01:3;
r = -1 + (1+1)*rand(1,201);
z1 = z1 + r;

x1 = x1';
y1 = y1';
z1 = z1';
label1 = ones(length(y1),1);

Tclust1 = [label1,x1,y1,z1];
%----------------------------
x2 = 4:0.01:6;
r = -1 + (1+1)*rand(1,201);
x2 = x2 + r;

y2 = 10:0.01:12;
r = -1 + (1+1)*rand(1,201);
y2 = y2 + r;

z2 = 10:0.01:12;
r = -1 + (1+1)*rand(1,201);
z2 = z2 + r;

x2 = x2';
y2 = y2';
z2 = z2';

label2 = 2.*ones(length(y2),1);

Tclust2 = [label2,x2,y2,z2];

%----------------------------
x3 = 5:0.01:7;
r = -1 + (1+1)*rand(1,201);
x3 = x3 + r;

y3 = 13:0.01:15;
r = -1 + (1+1)*rand(1,201);
y3 = y3 + r;

z3 = 13:0.01:15;
r = -1 + (1+1)*rand(1,201);
z3 = z3 + r;

x3 = x3';
y3 = y3';
z3 = z3';

label3 = 3.*ones(length(y3),1);

Tclust3 = [label3,x3,y3,z3];

%----------------------------
T = [Tclust1;Tclust2;Tclust3];

data = T(:,2:4);
labels = T(:,1);

%LOOCV
n = length(labels);
hpartition = cvpartition(n,'leaveout');

% figure()
output = zeros(n,3);

for qq = 1:n
    %training data
idx_Train = training(hpartition,qq);
data_Train = data(idx_Train,1:3);
labels_Train = labels(idx_Train);

    %test data
idx_Test = test(hpartition,qq);
data_Test = data(idx_Test,1:3);
labels_Test = labels(idx_Test);

%train the model on the training data
Mdlquad3D = fitcdiscr(data_Train,labels_Train,'Discrimtype','quadratic');

%apply to the testing data
pred_lbl = predict(Mdlquad3D,data_Test);

(qq/n)*100
end

iscorrect=output(:,1)==output(:,2);
accuracy_full = (sum(iscorrect)/length(iscorrect))*100;

%plot data with boundaries -----------------------------------------------
data = T(:,2:4);
labels = T(:,1);
xmin = min(T(:,2));
xmax = max(T(:,2));
ymin = min(T(:,3));
ymax = max(T(:,3));
zmin = min(T(:,4));
zmax = max(T(:,4));

%group 1 and group 2
K12 = Mdlquad3D.Coeffs(1,2).Const;
L12 = Mdlquad3D.Coeffs(1,2).Linear;
Q12 = Mdlquad3D.Coeffs(1,2).Quadratic;

%group 2 and group 3
K23 = Mdlquad3D.Coeffs(2,3).Const;
L23 = Mdlquad3D.Coeffs(2,3).Linear;
Q23 = Mdlquad3D.Coeffs(2,3).Quadratic;

%group 1 and group 3
K13 = Mdlquad3D.Coeffs(1,3).Const;
L13 = Mdlquad3D.Coeffs(1,3).Linear;
Q13 = Mdlquad3D.Coeffs(1,3).Quadratic;

figure()
scatter3(x1,y1,z1,'g')
hold on
scatter3(x2,y2,z2,'b')
scatter3(x3,y3,z3,'r')

%THESE ARE CORRECT FOR LDA - HOW TO INCLUDE Qxx? Is there a better way
%than fimplicit3?
f12 = @(x1,x2,x3) K12 + L12(1)*x1 + L12(2)*x2 + L12(3)*x3;
h2 = fimplicit3(f12,[xmin xmax ymin ymax zmin zmax]);
f23 = @(x1,x2,x3) K23 + L23(1)*x1 + L23(2)*x2 + L23(3)*x3;
h3 = fimplicit3(f23,[xmin xmax ymin ymax zmin zmax]);
f13 = @(x1,x2,x3) K13 + L13(1)*x1 + L13(2)*x2 + L13(3)*x3;
h4 = fimplicit3(f13,[xmin xmax ymin ymax zmin zmax]);

title('Full data set')
xlabel('x')
ylabel('y')
zlabel('z')
legend('Data 1','Data 2','Data 3','Boundary 1 -> 2','Boundary 2 -> 3','Boundary 1 -> 3','location','eastoutside') 
% 

供参考：

二维线性：

K = MdlLinear.Coeffs(1,2).Const;
L = MdlLinear.Coeffs(1,2).Linear;
f12 = @(x1,x2) K + L(1)*x1 + L(2)*x2;
h2 = fimplicit(f12,[xmin xmax ymin ymax]);

二维二次方：

K = MdlQuad.Coeffs(1,2).Const;
L = MdlQuad.Coeffs(1,2).Linear;
Q = MdlQuad.Coeffs(1,2).Quadratic;
f12 = @(x1,x2) K + L(1)*x1 + L(2)*x2 + Q(1,1)*x1.^2 + (Q(1,2)+Q(2,1))*x1.*x2 + Q(2,2)*x2.^2;
h2 = fimplicit(f12,[xmin xmax ymin ymax]);

Answer 1

从

classificationdiscriminant

文档中，我们有

Const

— 标量

Linear

— 具有 p 个分量的向量，其中 p 是 X
中的列数

Quadratic

— p×p 矩阵，存在于二次 DiscrimType

i类和j类的边界方程为

Const + Linear * x + x' * Quadratic * x = 0

其中 x 是长度为 p 的列向量。

您可以通过乘以并写出

x'Qx

来计算出二次项的所有矩阵系数，您已经为线性系数完成了（无论是否意识到）。

但是，将其视为矩阵方程并编写一个接受

x1,x2,x3

并使用矩阵评估边界的函数会更容易。

您可以通过创建一个本地函数

fKLQ(K,L,Q,x1,x2,x3)

来完成此操作，您可以将其与现有的匿名函数句柄

f12

等别名，如下所示：

f12 = @(x1,x2,x3) fKLQ( K12, L12, Q12, x1, x2, x3 );
h2 = fimplicit3(f12,[xmin xmax ymin ymax zmin zmax]);
f23 = @(x1,x2,x3) fKLQ( K23, L23, Q23, x1, x2, x3 );
h3 = fimplicit3(f23,[xmin xmax ymin ymax zmin zmax]);
f13 = @(x1,x2,x3) fKLQ( K13, L13, Q13, x1, x2, x3 );
h4 = fimplicit3(f13,[xmin xmax ymin ymax zmin zmax]);

function f = fKLQ( K, L, Q, x1, x2, x3 )
    x = [x1(:),x2(:),x3(:)];
    n = numel(x1);
    if n == 1
        x = x.';
    end
    f = zeros(1,n);
    for ii = 1:n
        f(ii) = K + L.'*x(:,ii) + x(:,ii).'*Q*x(:,ii);
    end
end

x1,x2,x3

的

fimplicit

输入要么是行向量，要么是标量，所以在正确形成

x

时必须处理它 - 可能有比我所展示的更简单的方法，但我相信它在至少这样工作。

输出将集群划分，如图所示。

注意：这里仍然有一些关于矢量化的警告，我相信这与我编写矩阵乘法的方式有关，MATLAB 警告评估输出可能不符合预期。您可以通过显式添加每个系数来避免这种情况，或者也许还有其他一些公式，但与此同时，我认为这是良性的。