数学建模常用代码

时间:2024-03-24 10:42:01

SVM分类器

1.命令函数部分:

clear;%清屏
clc;
X =load('data.txt');
n = length(X);%总样本数量
y = X(:,4);%类别标志
X = X(:,1:3);
TOL = 0.0001;%精度要求
C = 1;%参数,对损失函数的权重
b = 0;%初始设置截距b
Wold = 0;%未更新a时的W(a)
Wnew = 0;%更新a后的W(a)
for i = 1 : 50%设置类别标志为1或者-1
    y(i) = -1;
end
a = zeros(n,1);%参数a
for i = 1 : n%随机初始化a,a属于[0,C]
        a(i) = 0.2;
end

%为简化计算,减少重复计算进行的计算
K = ones(n,n);
for i = 1 :n%求出K矩阵,便于之后的计算
    for j = 1 : n
        K(i,j) = k(X(i,:),X(j,:));
    end
end
sum = zeros(n,1);%中间变量,便于之后的计算,sum(k)=sigma a(i)*y(i)*K(k,i);
for k = 1 : n
    for i = 1 : n
        sum(k) = sum(k) + a(i) * y(i) * K(i,k);
    end
end

while 1%迭代过程
    
%启发式选点
n1 = 1;%初始化,n1,n2代表选择的2个点
n2 = 2;
%n1按照第一个违反KKT条件的点选择
while n1 <= n
    if y(n1) * (sum(n1) + b) == 1 && a(n1) >= C && a(n1) <=  0
         break;
    end
    if y(n1) * (sum(n1) + b) > 1 && a(n1) ~=  0
           break;
    end
    if y(n1) * (sum(n1) + b) < 1 && a(n1) ~=C
          break;
    end
     n1 = n1 + 1;              
end
%n2按照最大化|E1-E2|的原则选取
E1 = 0;
E2 = 0;
maxDiff = 0;%假设的最大误差
E1 = sum(n1) + b - y(n1);%n1的误差
for i = 1 : n
    tempSum = sum(i) + b - y(i);
    if abs(E1 - tempSum)> maxDiff
        maxDiff = abs(E1 - tempSum);
        n2 = i;
        E2 = tempSum;
    end
end

%以下进行更新
a1old = a(n1);
a2old = a(n2);
KK = K(n1,n1) + K(n2,n2) - 2*K(n1,n2);
a2new = a2old + y(n2) *(E1 - E2) / KK;%计算新的a2
%a2必须满足约束条件
S = y(n1) * y(n2);
if S == -1
    U = max(0,a2old - a1old);
    V = min(C,C - a1old + a2old);
else
    U = max(0,a1old + a2old - C);
    V = min(C,a1old + a2old);
end
if a2new > V
    a2new = V;
end
if a2new < U
    a2new = U;
end
a1new = a1old + S * (a2old - a2new);%计算新的a1
a(n1) = a1new;%更新a
a(n2) = a2new;

%更新部分值
sum = zeros(n,1);
for k = 1 : n
    for i = 1 : n
        sum(k) = sum(k) + a(i) * y(i) * K(i,k);
    end
end
Wold = Wnew;
Wnew = 0;%更新a后的W(a)
tempSum = 0;%临时变量
for i = 1 : n
    for j = 1 : n
    tempSum= tempSum + y(i )*y(j)*a(i)*a(j)*K(i,j);
    end
    Wnew= Wnew+ a(i);
end
Wnew= Wnew - 0.5 * tempSum;
%以下更新b:通过找到某一个支持向量来计算
support = 1;%支持向量坐标初始化
while abs(a(support))< 1e-4 && support <= n
    support = support + 1;
end
b = 1 / y(support) - sum(support);
%判断停止条件
if abs(Wnew/ Wold - 1 ) <= TOL
    break;
end
end
%输出结果:包括原分类,辨别函数计算结果,svm分类结果
for i = 1 : n
    fprintf('第%d点:原标号 ',i);
    if i <= 50
        fprintf('-1');
    else
        fprintf(' 1');
    end
    fprintf('    判别函数值%f      分类结果',sum(i) + b);
    if abs(sum(i) + b - 1) < 0.5
        fprintf('1\n');
    else if abs(sum(i) + b + 1) < 0.5
            fprintf('-1\n');
        else
            fprintf('归类错误\n');
        end
    end
end

2.名为f的功能函数部分:
 

function y = k(x1,x2)
    y = exp(-0.5*norm(x1 - x2).^2);
end

K-means算法代码

function [Idx, Center] = K_means(X, xstart)
% K-means聚类
% Idx是数据点属于哪个类的标记,Center是每个类的中心位置
% X是全部二维数据点,xstart是类的初始中心位置

len = length(X);        %X中的数据点个数
Idx = zeros(len, 1);    %每个数据点的Id,即属于哪个类

C1 = xstart(1,:);       %第1类的中心位置
C2 = xstart(2,:);       %第2类的中心位置
C3 = xstart(3,:);       %第3类的中心位置

for i_for = 1:100
    %为避免循环运行时间过长,通常设置一个循环次数
    %或相邻两次聚类中心位置调整幅度小于某阈值则停止
    
    %更新数据点属于哪个类
    for i = 1:len
        x_temp = X(i,:);    %提取出单个数据点
        d1 = norm(x_temp - C1);    %与第1个类的距离
        d2 = norm(x_temp - C2);    %与第2个类的距离
        d3 = norm(x_temp - C3);    %与第3个类的距离
        d = [d1;d2;d3];
        [~, id] = min(d);   %离哪个类最近则属于那个类
        Idx(i) = id;
    end
    
    %更新类的中心位置
    L1 = X(Idx == 1,:);     %属于第1类的数据点
    L2 = X(Idx == 2,:);     %属于第2类的数据点
    L3 = X(Idx == 3,:);     %属于第3类的数据点
    C1 = mean(L1);      %更新第1类的中心位置
    C2 = mean(L2);      %更新第2类的中心位置
    C3 = mean(L3);      %更新第3类的中心位置
end

Center = [C1; C2; C3];  %类的中心位置


%演示数据
%% 1 random sample
%随机生成三组数据
a = rand(30,2) * 2;
b = rand(30,2) * 5;
c = rand(30,2) * 10;
figure(1);
subplot(2,2,1); 
plot(a(:,1), a(:,2), 'r.'); hold on
plot(b(:,1), b(:,2), 'g*');
plot(c(:,1), c(:,2), 'bx'); hold off
grid on;
title('raw data');

%% 2 K-means cluster
X = [a; b; c];  %需要聚类的数据点
xstart = [2 2; 5 5; 8 8];  %初始聚类中心
subplot(2,2,2);
plot(X(:,1), X(:,2), 'kx'); hold on
plot(xstart(:,1), xstart(:,2), 'r*'); hold off
grid on;
title('raw data center');

[Idx, Center] = K_means(X, xstart);
subplot(2,2,4);
plot(X(Idx==1,1), X(Idx==1,2), 'kx'); hold on
plot(X(Idx==2,1), X(Idx==2,2), 'gx');
plot(X(Idx==3,1), X(Idx==3,2), 'bx');
plot(Center(:,1), Center(:,2), 'r*'); hold off
grid on;
title('K-means cluster result');

disp('xstart = ');
disp(xstart);
disp('Center = ');
disp(Center);