MATLAB算法与模型学习笔记

最新推荐文章于 2023-07-06 21:21:27 发布

胆怯与勇敢

最新推荐文章于 2023-07-06 21:21:27 发布

阅读量970

点赞数 2

文章标签： matlab 算法学习

本文链接：https://blog.csdn.net/qq_52652401/article/details/129247939

版权

文章目录

一、主成分分析
二、整数规划
三、线性规划
四、图论
五、聚类分析
六、高斯拟合
七、非线性规划
八、多目标规划
九、典型相关分析
- 9.1 职业满意度典型相关分析
- 9.2 中国城市竞争力与基础设施的典型相关分析
十、插值与拟合

一、主成分分析

%% 主成分分析程序
% 主成分回归分析
clc, clear
load sn.txt
[m, n] = size(sn);
x0 = sn(:, 1 : n-1); y0 = sn(:, n);
hg1 = [ones(m, 1), x0] \ y0; % 计算普通最小二乘法回归系数
hg1 = hg1'; % 变成行向量显示回归系数,其中第一个分量是常数项
fprintf('y = %f', hg1(1)); % 显示普通最小二乘法回归结果
for i = 2 : n
    if hg1(i) > 0
        fprintf('+ %f * x%d', hg1(i), i-1); 
    else
        fprintf('%f * x%d', hg1(i), i-1);     
    end
end
fprintf('\n')
r = corrcoef(x0); % 计算相关系数矩阵
xd = zscore(x0); % 对设计矩阵进行标准化处理
yd = zscore(y0); % 对y0进行标准化处理
[vec1, lamda, rate] = pcacov(r); % vec1特征向量,lamda特征值,rate各主成分的贡献率
f = repmat(sign(sum(vec1)), size(vec1, 1), 1); % 构造与vec1同维度的元素为±1的矩阵
vec2 = vec1 .* f; % 修改特征向量的正负号,使得特征向量的所有和为正
contr = cumsum(rate); % 计算累计贡献率
df = xd * vec2; % 计算所有主成分的得分 
num = input('请选择主成分的个数：'); % 通过累计贡献率交互式选择主成分的个数
hg21 = df(:, 1 : num) \ yd; % 主成分变量的回归系数,这里由于数据标准化,回归方程的常数项为0
hg22 = vec2(:, 1 : num) * hg21; % 标准化回归方程的系数
% 计算原始变量回归方程的系数
hg23 = [mean(y0) - std(y0) * mean(x0) ./ std(x0) * hg22, std(y0) * hg22' ./ std(x0)];
fprintf('y = %f', hg23(1)); % 显示主成分回归结果
for i = 2 : n
    if hg23(i) > 0
        fprintf('+ %f * x%d', hg23(i), i-1); 
    else
        fprintf('%f * x%d', hg23(i), i-1);     
    end
end
fprintf('\n')
% 下面计算两种回归分析的剩余标准差
rmse1 = sqrt(sum((hg1(1) + x0 * hg1(2 : end)' - y0).^ 2) / (m -n)); % 拟合n个参数
rmse2 = sqrt(sum((hg23(1) + x0 * hg23(2 : end)' - y0).^ 2) / (m -num)); % 拟合num个参数

二、整数规划

2.1 入门题

clear; clc
%% 第一题
C = [-20, -10]';
intcon = [1, 2]; % 限制第一个变量和第二个变量为整数
A = [5, 4;
     2, 5];
b = [24, 13]';
lb = [0, 0]';
[x1, fval1] = intlinprog(C, intcon, A, b, [], [], lb); 
fval1 = -fval1; % 最大化问题

%% 第二题
C = [18, 23, 5]';
intcon = 3; % 限制第三个变量为整数
A = [107, 500, 0;
     72, 121, 65;
     -107, -500, 0;
     -72, -121, -65];
b = [50000, 2250, -500, -2000]';
lb = [0, 0, 0]';
[x2, fval2] = intlinprog(C, intcon, A, b, [], [], lb);

%% 第三题(0-1规划)
C = [-3, -2, -1]';
intcon = 3; % 表示第三个变量为整数
A = [1, 1, 1];
b = 7;
Aeq = [4, 2, 1];
beq = 12;
lb = [0, 0, 0]';
ub = [+inf, +inf, 1]'; % 表示0-1规划
[x3, fval3] = intlinprog(C, intcon, A, b, Aeq, beq, lb, ub);

2.2 典型例题

2.2.1 背包问题(0-1规划)

clear; clc
%% 背包问题(0-1规划)
C = [540, 200, 180, 350, 60, 150, 280, 450, 320, 120]';
intcon = 1 : 10;
A = [6, 3, 4, 5, 1, 2, 3, 5, 4, 2];
b = 30;
lb = zeros(10, 1);
ub = ones(10, 1);
[x1, fval1] = intlinprog(-C, intcon, A, b, [], [], lb, ub); 
fval1 = -fval1;

2.2.2 指派问题(0-1规划)

%% 指派问题(0-1规划)
t = [66.8, 75.6, 87, 58.6;
     57.2, 66, 66.4, 53;
     78, 67.8, 84.6, 59.4;
     70, 74.2, 69.6, 57.2;
     67.4, 71, 83.8, 62.4];
C = zeros(20, 1);
for i = 1 : 5
    for j = 1 : 4
        C((i - 1) * 4 + j) = t(i, j); 
    end
end
intcon = 1 : 20;
A = zeros(5, 20);
for i = 1 : 5
   for j = 1 : 4
       A(i, (i - 1) * 4 + j) = 1;
   end
end
b = ones(5, 1);
Aeq = zeros(4, 20);
% Aeq = [repmat(eye(4), 1, 5)]; % 简洁
for i = 1 : 4
   for j = 1 : 5
       Aeq(i, i + 4 * (j - 1)) = 1;
   end
end
beq = ones(4, 1);
lb = zeros(20, 1);
ub = ones(20, 1);
[x2, fval2] = intlinprog(C, intcon, A, b, Aeq, beq, lb, ub); 
x2 = reshape(x2, 4, 5)';

2.2.3 钢管切割问题(混合使用切割方案)

%% 钢管切割问题(混合使用切割方案)
% 找出所有的切割方案(枚举法)
for i = 0 : 2
    for j = 0 : 3
        for k = 0 : 6
            len = 2.9 * i + 2.1 * j + k;
            if len <= 6.9 && len >= 6
                disp([i, j, k])
            end
        end
    end
end
C = ones(7, 1);
intcon = 1 : 7;
A = [-1, -2, 0, 0, 0, 0, -1;
     0, 0, -3, -2, -1, 0, -1;
     -4, -1, 0, -2, -4, -6, -1];
b = [-100, -100, -100]';
lb = zeros(7, 1);
[x3, fval3] = intlinprog(C, intcon, A, b, [], [], lb);

三、线性规划

3.1 入门题

clear; clc
%% 第一题
C = [-5, -4, -6]'; % 与标准型保持一致,不转置也行
A = [1, -1, 1;
     3, 2, 4;
     3, 2, 0];
b = [20, 42, 30]';
lb = [0, 0, 0]';
[x1, fval1] = linprog(C, A, b, [], [], lb); % 没有ub表示没有上界约束

%% 第二题(此题有多解)
C = [0.04, 0.15, 0.1, 0.125]'; % 与标准型保持一致,不转置也行
A = [-0.03, -0.3, 0, -0.15;
     0.14, 0, 0, 0.07];
b = [-32, 42]';
Aeq = [0.05, 0, 0.2, 0.1];
beq = 24;
lb = [0, 0, 0, 0]';
[x2, fval2] = linprog(C, A, b, Aeq, beq, lb); % 没有ub表示没有上界约束

%% 第三题
C = [-2, -3, 5]'; % 与标准型保持一致,不转置也行
A = [-2, 5, -1;
     1, 3, 1];
b = [-10, 12]';
Aeq = [1, 1, 1];
beq = 7;
lb = [0, 0, 0]';
[x3, fval3] = linprog(C, A, b, Aeq, beq, lb); % 没有ub表示没有上界约束
fval3 = -fval3; % 结果取最大值,需添加符号

%% 多个解的情况
% 例如:min z = x1 + x2  s.t. x1 + x2 >= 10
C = [1, 1]';
A = [-1, -1];
b = -10;
[x4, fval4] = linprog(C, A, b); % 有多个解时,返回一个值

%% 不存在解的情况
% 例如:min z = x1 + x2  s.t. x1 + x2 = 10、x1 + 2 * x2 <=8
C = [1, 1]';
A = [1, 2];
b = 8;
Aeq = [1, 1];
beq = 10;
lb = [0, 0]';
[x5, fval5] = linprog(C, A, b, Aeq, beq, lb);

3.2 典型例题

3.2.1 生产决策问题

clear; clc
%% 生产决策问题(本题应为整数规划,仅供入门参考)
format long g % 计算结果显示为长数字
% 目标函数系数向量
C = zeros(9, 1);
C(1) = 0.75; C(2) = 0.7753; C(3) = -0.375;
C(4) = -783 / 1750; C(5) = -0.35; C(6) = -0.5;
C(7) = -0.2889; C(8) = 1.15; C(9) = 1.9148 - 8613 / 7000;
% 系数矩阵(A) 不等式约束
A = zeros(5, 9);
A(1, 1) = 1; A(1, 6) = 2;
A(2, 2) = 7; A(2, 7) = 9; A(2, 9) = 12;
A(3, 3) = 3; A(3, 8) = 4;
A(4, 4) = 4; A(4, 9) = 11;
A(5, 5) = 7; 
% 常数项(b) 不等式约束
b = zeros(5, 1);
b(1) = 1200;
b(2) = 10000;
b(3) = 2000;
b(4) = 7000;
b(5) = 4000;
% 系数矩阵(Aeq) 等式约束
Aeq  = zeros(2, 9);
Aeq(1, 1) = 1; Aeq(1, 2) = 1; Aeq(1, 3) = -1; Aeq(1, 4) = -1; Aeq(1, 5) = -1;
Aeq(2, 6) = 1; Aeq(2, 7) = 1; Aeq(2, 8) = -1; 
% 常数项(beq) 等式约束
beq = zeros(2, 1);
% 下界(lb)
lb = zeros(9, 1);
[x1, fval1] = linprog(-C, A, b, Aeq, beq, lb);
fval1 = -fval1; % 题目为求最大值

3.2.2 生成决策问题（整数规划）

% 将其改为整数规划
intcon = 1 : 9;
[x, fval] = intlinprog(-C, intcon, A, b, Aeq, beq, lb);
fval = -fval;
%% 投料问题(运输问题)
format long g % 计算结果显示为长数字
% 目标函数系数向量
M = [5, 1;
     2, 7]; % 产地
N = [1.25, 1.25;
     8.75, 0.75;
     0.5, 4.75;
     5.75, 5;
     3, 6.5;
     7.25, 7.25]; % 销地
C = zeros(12, 1);
for i = 1 : length(M)
   for j = 1 : length(N)
      distance = (M(i, 1) - N(j, 1)) ^ 2 + (M(i, 2) - N(j, 2)) ^ 2;
      C((i - 1) * length(N) + j) = sqrt(distance);
   end
end
% 系数矩阵(A) 不等式约束
A = zeros(2, 12);
A(1, 1 : 6) = 1; % 简洁
A(2, 7 : 12) = 1; % 简洁
% A(1, 1) = 1; A(1, 2) = 1; A(1, 3) = 1; A(1, 4) = 1; A(1, 5) = 1; A(1, 6) = 1;
% A(2, 7) = 1; A(2, 8) = 1; A(2, 9) = 1; A(2, 10) = 1; A(2, 11) = 1; A(2, 12) = 1;
% 常数项(b) 不等式约束
b = [20, 20]';
% 系数矩阵(Aeq) 等式约束
Aeq  = zeros(6, 12);
for i = 1 : 6
   Aeq(i, i) = 1; Aeq(i, i + 6) = 1; 
end
Aeq = [eye(6), eye(6)]; % 简洁
% Aeq(1, 1) = 1; Aeq(1, 7) = 1; 
% Aeq(2, 2) = 1; Aeq(2, 8) = 1; 
% Aeq(3, 3) = 1; Aeq(3, 9) = 1;
% Aeq(4, 4) = 1; Aeq(4, 10) = 1; 
% Aeq(5, 5) = 1; Aeq(5, 11) = 1; 
% Aeq(6, 6) = 1; Aeq(6, 12) = 1;
% 常数项(beq) 等式约束
beq = [3, 5, 4, 7, 6, 11]';
% 下界(lb)
lb = zeros(12, 1);
[x2, fval2] = linprog(C, A, b, Aeq, beq, lb);
x2  = reshape(x2, 6, 2)'; % 展示效果更加美观

四、图论

4.1 最短路径

%% 最短路径
% Dijkstra算法
shortestpath();
% Floyd算法:求图中任意两点之间最短路径
distances()

4.2 最小生成树

4.2.1 prim算法

%% 最小生成树
% prim算法
clc; clear
a = zeros(7);
a(1, 2) = 50; a(1, 3) = 60; a(2, 4) = 65;
a(2, 5) = 40; a(3, 4) = 52; a(3, 7) = 45;
a(4, 5) = 50; a(4, 6) = 30; a(4, 7) = 42;
a(5, 6) = 70;
a = a + a'; a(a == 0) = inf;
result = []; p = 1; tb = 2 : length(a);
while size(result, 2) ~= length(a) - 1
    temp = a(p, tb); temp = temp(:);
    d = min(temp);
    [jb, kb] = find(a(p, tb) == d, 1); % 找第一个最小的值
    j = p(jb); k = tb(kb);
    result = [result, [j; k; d]];
    p = [p, k];
    tb(tb == k) = [];
end
result

4.2.2 Kruskal算法

% Kruskal算法
clc; clear
% 这里给出邻接矩阵的另外一种输入方式
a(1, [2, 3]) = [50, 60]; a(2, [4, 5]) = [65, 40];
a(3, [4, 7]) = [52, 45]; a(4, [5, 6]) = [50, 30];
a(4, 7) = 42; a(5, 6) = 70;
[i, j, b] = find(a);
data = [i'; j'; b']; index = data(1 : 2, :);
loop = length(a) - 1;
result = [];
while length(result) < loop
    temp = min(data(3, :));
    flag = find(data(3, :) == temp);
    flag = flag(1);
    v1 = index(1, flag); v2 = index(2, flag);
    if v1 ~= v2
        result = [result, data(:, flag)];
    end
    index(index == v2) = v1;
    data(:, flag) = [];
    index(:, flag) = [];
end
result

4.3 Matlab的图论工具箱

%% Matlab的图论工具箱
% 无向图最短路径及长度
clc; clear
a(1, 2) = 2; a(1, 3) = 8; a(1, 4) = 1; a(2, 3) = 1;
a(2, 3) = 6; a(2, 5) = 1; a(3, 4) = 7; a(3, 5) = 5;
a(3, 6) = 1; a(3, 7) = 2; a(4, 7) = 9; a(5, 6) = 3;
a(5, 8) = 2; a(5, 9) = 9; a(6, 7) = 4; a(6, 8) = 6;
a(7, 9) = 3; a(7, 10) = 1; a(8, 9) = 7; a(8, 11) = 9;
a(9, 10) = 1; a(9, 11) = 2; a(10, 11) = 4;
a = a'; % Matlab工具箱要求数据是下三角矩阵
[i, j, v] = find(a);
b = sparse(i, j, v, 11, 11); % 构造系数矩阵
% Directed是标志图为有向或无向的属性
[x, y, z] = graphshortestpath(b, 1, 11, 'Directed', false)

4.4 渡河问题

% 渡河问题
clc; clear
a = [1, 1, 1, 1;
     1, 1, 1, 0;
     1, 1, 0, 1;
     1, 0, 1, 1;
     1, 0, 1, 0;
     0, 1, 0, 1;
     0, 1, 0, 0;
     0, 0, 1, 0;
     0, 0, 0, 1;
     0, 0, 0, 0]; % 每一行是一个可行状态
b = [1, 0, 0, 0;
     1, 1, 0, 0;
     1, 0, 1, 0;
     1, 0, 0, 1]; % 每一行是一个可以转移状态
w = zeros(10); % 邻接矩阵初始化
for i = 1 : 9
    for j = i + 1 : 10
        for k = 1 : 4
            if strfind(xor(a(i, :), b(k, :)), a(j, :))
                w(i,j) = 1; 
            end 
        end
    end
end
w = w';  % 变成下三角形矩阵
c = sparse(w); % 构造稀疏矩阵
% 该图是无向图,Directed属性值为0
[x, y, z] = graphshortestpath(c, 1, 10, 'Directed', 0)
% 画出无向图
h = view(biograph(c, [], 'ShowArrows', 'off', 'ShowWeights', 'off'));
Edges = getedgesbynodeid(h); % 提取句柄h中的边集
set(Edges, 'LineColor', [0, 0, 0]); % 为了将来打印清楚,边画成黑色
set(Edges, 'LineWidth', 1.5); % 线型宽度设置为1.5

4.5 有向图最短路径及长度

% 有向图最短路径及长度
clc; clear
a = zeros(7);
a(1, 2) = 4; a(1, 3) = 2; a(2, 3) = 3;
a(2, 4) = 2; a(2, 5) = 6; a(3, 4) = 5;
a(3, 6) = 4; a(4, 5) = 2; a(4, 6) = 7;
a(5, 6) = 4; a(5, 7) = 8; a(6, 7) = 3;
b = sparse(a); % 构造稀疏矩阵
% 有向图Directed属性值为真或1,方法属性的默认值是Dijkstra
[x, y, z] = graphshortestpath(b, 1, 7, 'Directed', 1, 'Method', 'Bellman-Ford')
view(biograph(b, []))

4.6 最小生成树（工具箱）

% 最小生成树
clc; clear
x = [0, 5, 16, 20, 33, 23, 35, 25, 10];
y = [15, 20, 24, 20, 25, 11, 7, 0, 3];
xy = [x; y];
d = mandist(xy); % 求xy的两两列向量间的绝对值距离
d = tril(d); % 截取Matlab工具箱要求的下三角形矩阵
b = sparse(d); % 转化为稀疏矩阵
[ST, pred] = graphminspantree(b, 'Method', 'Kruskal') % 调用最小生成树的命令
st = full(ST); % 把最小生成树的稀疏矩阵转化为普通矩阵
TreeLength = sum(sum(st)) % 求最小生成树长度
view(biograph(ST, [], 'ShowArrows', 'off')) % 画出最小生成树

4.7 最大流

% 最大流
% 只能解决权重都为正值,且两个顶点之间不能有两条弧
% 添加虚拟顶点9,解决顶点3、4之间的两条弧
clc; clear
a = zeros(9);
a(1, 2) = 6; a(1, 3) = 4; a(1, 4) = 5; a(2, 3) = 3;
a(2, 5) = 9; a(2, 6) = 9; a(3, 4) = 4; a(3, 5) = 6;
a(3, 6) = 7; a(3, 7) = 3; a(4, 7) = 5; a(4, 9) = 2;
a(5, 8) = 12; a(6, 5) = 8; a(6, 8) = 10; a(7, 6) = 4;
a(7, 8) = 15; a(9, 3) = 2;
b = sparse(a);
[x, y, z] = graphmaxflow(b, 1, 8)
view(biograph(b, []))

4.8 旅行商问题

%% 旅行商问题
% 修改圈近似算法
clc; clear
a = zeros(6);
a(1, 2) = 56; a(1, 3) = 35; a(1, 4) = 21; a(1, 5) = 51;
a(1, 6) = 60; a(2, 3) = 21; a(2, 4) = 57; a(2, 5) = 78;
a(2, 6) = 70; a(3, 4) = 36; a(3, 5) = 68; a(3, 6) = 68;
a(4, 5) = 51; a(4, 6) = 61; a(5, 6) = 13;
a = a + a'; L = size(a, 1);
c = [5, 1 : 4, 6, 5]; % 选取初始圈
for k = 1 : L
    flag = 0; % 退出标志
    for m = 1 : L - 2 % m为算法中的i
        for n = m + 2 : L % n为算法中的j
            if a(c(m), c(n)) + a(c(m + 1), c(n + 1)) < a(c(m), c(m + 1)) + a(c(n), c(n + 1))
                c(m + 1 : n) = c(n : -1 : m + 1); flag = flag + 1; % 修改一次,标志加1
            end
        end
    end
    if flag == 0 % 一条边也没有修改,就返回
        long = 0; % 圈长的初始值
        for i = 1 : L
            long = long + a(c(i), c(i + 1)); % 求改良圈的长度
        end
        circle = c; % 返回修改圈
        return
    end
end

五、聚类分析

5.1 DBSCAN算法

%% DBSCAN算法
clc; clear;
load mydata; % 加载数据
epsilon = 0.5; % 半径
MinPts = 10; % 点数
C = 0; % 记录分的类数
n = size(X, 1); % X为加载的数据集
IDX = zeros(n, 1);  % 初始化全部为0,即全部为噪音点
D = pdist2(X, X); % 计算每个点之间的距离
visited = false(n, 1); % 是否已经访问
isnoise = false(n, 1); % 是否为噪音点
for i = 1 : n
    if ~visited(i)
        visited(i) = true; % 更新为已访问
        % 记录该访问点的临近点
        Neighbors = find(D(i, :) <= epsilon);
        if numel(Neighbors) < MinPts % numel()计算数组中元素数量
            % 附近点的数量小于已给数量,记为噪声点
            isnoise(i) = true;
        else
            C = C + 1;
            IDX(i) = C;
            % 开始搜索其临近点
            k = 1;
            while true
                j = Neighbors(k);
                if ~visited(j)
                    visited(j) = true;
                    Neighbors2 = find(D(j,:) <= epsilon);
                    if numel(Neighbors2) >= MinPts
                        Neighbors=[Neighbors Neighbors2];
                    end
                end
                if IDX(j) == 0
                    IDX(j) = C;
                end
                k = k + 1;
                if k > numel(Neighbors)
                    break;
                end
            end
        end
    end
end
% 如果只有两个指标的,就可以绘制图形
k = max(IDX);
Colors = hsv(k);
Legends = {};
for i = 0 : k
    Xi = X(IDX == i, :); % 第i类的所有值
    if i ~= 0
        Style = 'x';
        MarkerSize = 8;
        Color = Colors(i,:);
        Legends{end + 1} = ['Cluster' num2str(i)];
    else % 噪声点
        Style = 'o';
        MarkerSize = 6;
        Color = [0 0 0];
        if ~isempty(Xi)
            Legends{end + 1} = 'Noise';
        end
    end
    if ~isempty(Xi)
        plot(Xi(:,1), Xi(:,2), Style, 'MarkerSize', MarkerSize, 'Color', Color);
    end
    hold on;
end
hold off;
axis equal;
grid on; % 增加网格线
legend(Legends); % 标签

5.2 聚类分析程序

%% 聚类分析程序
clc, clear
a = [1, 0; 1, 1; 3, 2; 4, 3; 2, 5];
[m, n] = size(a);
d = mandist(a'); % mandist求矩阵列向量组之间的两两绝对值距离
d = tril(d); % 截取下三角形
nd = nonzeros(d); % 去掉d中的零元素,非零元素按列排列
d = union([], nd); % 去掉重复的非零元素
for i = 1 : m - 1
    nd_min = min(nd);
    [row, col] = find(d == nd_min);
    tm = union(row, col); % row和col归为一类
    tm = reshape(tm, 1, length(tm)); % 把数组tm变成行向量
    fprintf('第%d次合成,平台高度为%d时的分类结果为：%s\n', i, nd_min, int2str(tm));
    nd(nd == nd_min) = []; % 删除已经归类的元素
    if length(nd) == 1
        break 
    end
end

5.3 聚类分析（matlab统计工具箱）

% 使用matlab统计工具箱
clc, clear
a = [1, 0; 1, 1; 3, 2; 4, 3; 2, 5];
y = pdist(a, 'cityblock'); % 求a的两两行向量间的绝对值距离
yc = squareform(y); % 变换成距离矩阵
z = linkage(y); % 产生等级聚类书
dendrogram(z) % 画聚类图
T = cluster(z, 'maxclust', 3); % 把对象划分成3类
% zsore(X) 对数据矩阵进行标准化处理
for i = 1 : 3
    tm = find(T == i); % 求第i类的对象
    tm = reshape(tm, 1, length(tm)); % 变成行向量
    fprintf('第%d类的有%s\n', i, int2str(tm)); %显示分类结果
end

5.4 变量聚类

% 对变量聚类
clc, clear
a = textread('ch.txt');
d = 1 - abs(a); % 进行数据变换,把相关系数转化为距离
d = tril(d); % 提取出d矩阵的下三角部分
b = nonzeros(d); % 去掉d中的零元素
b = b'; % 化为行向量(对变量聚类)
z = linkage(b, 'complete'); % 按最长距离法聚类
y = cluster(z, 'maxclust', 2); % 把变量划分为两类
ind1 = find(y == 1); ind1 = ind1'; % 显示第一类对应的变量标号
ind2 = find(y == 2); ind2 = ind2'; % 显示第二类对应的变量标号
h = dendrogram(z); % 画聚类图
set(h, 'Color', 'k', 'LineWidth', 1.3) % 把聚类图线的颜色改为黑色,线宽加粗

5.5 我国各地区普通高等教育发展情况分析

5.5.1 R型聚类(选取指标)

    % R型聚类(选取指标)
clc, clear
a = load('gj.txt'); % 把原始数据保存在纯文本文件gj.txt中
b = zscore(a); % 数据标准化
r = corrcoef(b); % 计算相关系数矩阵
% d = tril(1 - r); d = nonzeros(d)'; % 另一种计算距离方法
d = pdist(b', 'correlation'); % 计算相关系数到处的距离
z = linkage(d, 'average'); % 按类平均法聚类
h = dendrogram(z); % 画聚类图
set(h, 'Color', 'k', 'LineWidth', 1.3) % 把聚类图线的颜色改为黑色,线宽加粗
T = cluster(z, 'maxclust', 6); % 把变量划分为6类
for i = 1 : 6
    tm = find(T == i); % 求第i类对象
    tm = reshape(tm, 1, length(tm)); % 变成行向量
    fprintf('第%d类的有%s\n', i, int2str(tm)); % 显示分类结果
end

5.5.2 Q型聚类

    % Q型聚类
clc, clear
load gj.txt  % 把原始数据保存在纯文本文件gj.txt中
gj(:, 3 : 6) = []; % 删除数据矩阵的第3-6列,即使用变量1,2,7,8,9,10
gj = zscore(gj); % 数据标准化
y = pdist(gj); % 求对象间的欧氏距离,每行是一个对象
z = linkage(y, 'average'); % 按类平均法聚类
h = dendrogram(z); % 画聚类图
set(h, 'Color', 'k', 'LineWidth', 1.3) % 把聚类图线的颜色改为黑色,线宽加粗 
for k = 3 : 5
    fprintf('划分为%d类的结果如下：\n', k)
    T = cluster(z, 'maxclust', k); % 把样本点划分为k类
    for i = 1 : k
        tm = find(T == i); % 求第i类对象
        tm = reshape(tm, 1, length(tm)); % 变成行向量
        fprintf('第%d类的有%s\n', i, int2str(tm)); % 显示分类结果
    end
    if k == 5
        break
    end
    fprintf('* * * * * * * * * * * * * * * * * * * * * * * * * * * *\n')
end

六、高斯拟合

clc; clear
y = [108.68, 118.71, 120.13, 106.71, 88.5, 95.15, 98.85, 89.11, 74.33, 77.29]';
x = (1 : 2 : 2 * length(y))';
%% fitgp matlab自带的高斯拟合
gpr = fitrgp(x, y, 'Sigma', 0.1);
x_test = linspace(1, 20, 100)';
[y_test, ~ , limit] = predict(gpr, x_test);
%% 画图
h1 = fill([x_test', fliplr(x_test')], [limit(:, 1)', fliplr(limit(:, 2)')], 'r', 'DisplayName', 'uncertain');
hold on
COLOR = [184  184  246] / 255;
COLOR1 = [207  155  255] / 255;
COLOR2 = [190  210  254] / 255;
COLOR3 = [184  184  246] / 255;
h1.FaceColor = COLOR; % 定义区间的填充颜色      
h1.EdgeColor = [1, 1, 1]; % 边界颜色设置为白色
alpha .2  % 设置透明色
plot(x_test, y_test, 'Color', COLOR1, 'LineWidth', 1, 'DisplayName', 'prediction') 
plot(x_test, limit(:, 1), '--', 'Color', COLOR2, 'DisplayName', 'Lower  Limit')
plot(x_test, limit(:, 2), '--', 'Color', COLOR3, 'DisplayName', 'Upper  Limit')
plot(x, y, '+', 'DisplayName', 'True Data')
legend('show', 'Location', 'Best')

七、非线性规划

7.1 入门题

%% 非线性规划
clear; clc
%% 第一问
format long g
x0 = [0, 0]; % 给定的任意初始值
A = [-2, 3]; b = 6; % 可以将线性约束改为非线性约束,但不推荐使用
[x1, fval1] = fmincon(@fun1, x0, A, b, [], [], [], [], @nonlfun1);
fval1 = -fval1; % 求极大值
%% 其他算法(初始值选取非常重要)
% 使用interior point算法(内点法)
option = optimoptions('fmincon', 'Algorithm', 'interior-point');
[x2, fval2] = fmincon(@fun1, x0, A, b, [], [], [], [], @nonlfun1, option);
fval2 = -fval2;
% 使用SQP算法(序列二次规划法)
option = optimoptions('fmincon', 'Algorithm', 'sqp');
[x3, fval3] = fmincon(@fun1, x0, A, b, [], [], [], [], @nonlfun1, option);
fval3 = -fval3; % 结果较差
x00 = [1, 1]; % 更改初始值
[x4, fval4] = fmincon(@fun1, x00, A, b, [], [], [], [], @nonlfun1, option);
fval4 = -fval4; % 效果变好,说明序列二次规划法适应性较差
% 使用active set算法(有效集法)
option = optimoptions('fmincon', 'Algorithm', 'active-set');
[x5, fval5] = fmincon(@fun1, x0, A, b, [], [], [], [], @nonlfun1, option);
fval5 = -fval5;
% 使用trust-region-reflective算法(信赖域反射算法)
option = optimoptions('fmincon', 'Algorithm', 'trust-region-reflective');
[x6, fval6] = fmincon(@fun1, x0, A, b, [], [], [], [], @nonlfun1, option);
fval6 = -fval6;

%% 第二问
format long g
x0 = [0, 0, 0];
lb = [0, 0, 0]';
[x7, fval7] = fmincon(@fun2, x0, [], [], [], [], lb, [], @nonlfun2);

%% 第三问
format long g
x0 = [22.58, 12.58, 12.13]; % 可以使用蒙特卡罗模拟得到初始解
A = [1, -2, -2;
     1, 2, 2];
b = [0, 72]';
Aeq = [1 ,-1, 0];
beq = 10;
lb = [-inf, 10, -inf]';
ub = [+inf, 20, +inf]';
[x8, fval8] = fmincon(@fun3, x0, A, b, Aeq, beq, lb, ub);
fval8 = -fval8; % 求极大值

7.2 典型题

7.2.1 选址问题

%% 典型例题
clear; clc
%% 选址问题
format long g
% 系数矩阵(A) 不等式约束
A = zeros(2, 16);
A(1, 1 : 6) = 1; % 简洁
A(2, 7 : 12) = 1; % 简洁
% 常数项(b) 不等式约束
b = [20, 20]';
% 系数矩阵(Aeq) 等式约束
Aeq = [eye(6), eye(6), zeros(6, 4)]; % 简洁
% 常数项(beq) 等式约束
beq = [3, 5, 4, 7, 6, 11]';
% 下界(lb)
lb = zeros(16, 1);
% 进行求解
x0 = [3, 5, 0, 7, 0, 1, 0, 0, 4, 0, 6, 10, 5, 2, 1, 7]; %利用第一问的答案求解 
[x1, fval1] = fmincon(@fun4, x0, A, b, Aeq,beq, lb);
x1 = reshape(x1(1 : 12), 6, 2)';

7.2.2 飞行管理问题

%% 飞行管理问题
format long g
% 画六架飞机的位置
figure(1) % 生成一个图形
box on % 显示为封闭的盒子
% 绘制飞机的初始位置
data = [150, 140, 243;
        85, 85, 236;
        150, 155, 220.5;
        145, 50, 159;
        130, 150, 230;
        0, 0, 52]; % x坐标 y坐标 方向角
plot(data(:, 1), data(:, 2), '.r')
axis([0, 160, 0, 160]) % 设置坐标轴刻度范围
hold on % 接着绘制
% 在图上标上注释
for i = 1 : 6
    txt = ['飞机', num2str(i)];
    text(data(i, 1) + 2, data(i, 2) + 2, txt, 'FontSize', 8)
end
% 求解非线性规划问题
x0 = [0, 0, 0, 0, 0, 0]; % 飞机调整角度的初始值
lb = -pi / 6 * ones(6, 1); % 下界
ub = pi / 6 * ones(6, 1); % 上界
[x2, fval2] = fmincon(@fun5, x0, [], [], [], [], lb, ub, @nonlfun5);
x2 = x2 * pi / 180; % 将弧度制转换为角度制

八、多目标规划

%% 第一种方式
a = [-1, -1, 0, 0;
     0, 0, -1, -1;
     3, 0, 2, 0;
     0, 3, 0, 2];
b = [-30, -30, 120, 48]';
c1 = [-100, -90, -80, -70]';
c2 = [0, 3, 0, 2]';
% 第一个目标函数的目标值
[x1, g1] = linprog(c1, a, b, [], [], zeros(4, 1));
% 第二个目标函数的目标值
[x2, g2] = linprog(c2, a, b, [], [], zeros(4, 1));
g3 = [g1, g2]'; % 目标goal的值
% 这里的权重weight=目标goal的绝对值
[x, fval] = fgoalattain('fun', rand(4, 1), g3, abs(g3), a, b, [], [], zeros(4, 1));

%% 第二种方式
clc; clear
a = [-1, -1, 0, 0;
     0, 0, -1, -1;
     3, 0, 2, 0;
     0, 3, 0, 2];
b = [-30, -30, 120, 48]';
c1 = [-100, -90, -80, -70];
c2 = [0, 3, 0, 2];
Fun = @(x)[c1; c2] * x; % 用匿名函数定义目标向量
% 第一个目标函数的目标值
[x1, g1] = linprog(c1, a, b, [], [], zeros(4, 1));
% 第二个目标函数的目标值
[x2, g2] = linprog(c2, a, b, [], [], zeros(4, 1));
g3 = [g1, g2]'; % 目标goal的值
% 这里的权重weight=目标goal的绝对值
[x, fval] = fgoalattain(Fun, rand(4, 1), g3, abs(g3), a, b, [], [], zeros(4, 1));

九、典型相关分析

9.1 职业满意度典型相关分析

%% 典型相关分析
% 职业满意度典型相关分析
clc; clear
load r.txt
n1 = 5; n2 = 7; num = min(n1, n2);
s11 = r(1 : n1, 1 : n1); % 提出X与X的相关系数
s12 = r(1 : n1, n1 + 1 : end); % 提出X与Y的相关系数
s21 = s12'; % 提出Y与X的相关系数
s22 = r(n1 + 1 : end, n1 + 1 : end); % 提出Y与Y的相关系数
m1 = inv(s11) * s12 * inv(s22) * s21; % 计算矩阵M1
m2 = inv(s22) * s21 * inv(s11) * s12; % 计算矩阵M2
[vec1, val1] = eig(m1); % 求M1的特征向量和特征值
for i = 1 : n1
    vec1(:, i) = vec1(:, i) / sqrt(vec1(:, i)' * s11 * vec1(:, i)); % 特征向量归一化,满足a's11a = 1
    vec1(:, i) = vec1(:, i) / sign(sum(vec1(:, i))); % 保证所有分量和为正
end
val1 = sqrt(diag(val1)); % 计算特征值的平方根
[val1, ind1] = sort(val1, 'descend'); % 按照从大到小排列
a = vec1(:, ind1(1 : num)); % 取出X组的系数阵
dcoef1 = val1(1 : num); % 提出典型相关系数
xlswrite('bk.xls', a, 'Sheet1', 'A1') % 把计算结果写到Excel文件中去
flag = n1 + 2; % % 把计算结果写到Excel中的行计数变量
str = char(['A', int2str(flag)]); % str为Excel中写数据的起始位置
xlswrite('bk.xls', dcoef1', 'Sheet1', str)
[vec2, val2] = eig(m2); % 求M2的特征向量和特征值
for i = 1 : n2
    vec2(:, i) = vec2(:, i) / sqrt(vec2(:, i)' * s22 * vec2(:, i)); % 特征向量归一化,满足a's22a = 1
    vec2(:, i) = vec2(:, i) / sign(sum(vec2(:, i))); % 保证所有分量和为正
end
val2 = sqrt(diag(val2)); % 计算特征值的平方根
[val2, ind2] = sort(val2, 'descend'); % 按照从大到小排列
b = vec2(:, ind2(1 : num)); % 取出Y组的系数阵
dcoef2 = val2(1 : num); % 提出典型相关系数
flag = flag + 2; str = char(['A', int2str(flag)]); % str为Excel中写数据的起始位置
xlswrite('bk.xls', b, 'Sheet1', str) % 把计算结果写到Excel文件中去
flag = flag + n2 + 1; str = char(['A', int2str(flag)]); % str为Excel中写数据的起始位置
xlswrite('bk.xls', dcoef2', 'Sheet1', str)
x_u_r = s11 * a; % x,u的相关系数
y_v_r = s22 * b; % y,v的相关系数
x_v_r = s12 * b; % x,v的相关系数
y_u_r = s21 * a; % y,u的相关系数
flag = flag + 2; str = char(['A', int2str(flag)]);
xlswrite('bk.xls', x_u_r', 'Sheet1', str)
flag = flag + n1 + 1; str = char(['A', int2str(flag)]);
xlswrite('bk.xls', y_v_r', 'Sheet1', str)
flag = flag + n2 + 1; str = char(['A', int2str(flag)]);
xlswrite('bk.xls', x_v_r', 'Sheet1', str)
flag = flag + n1 + 1; str = char(['A', int2str(flag)]);
xlswrite('bk.xls', y_u_r', 'Sheet1', str)
mu = sum(x_u_r .^ 2) / n1 % X组原始变量被u_i解释的方差比例
mv = sum(x_v_r .^ 2) / n1 % X组原始变量被v_i解释的方差比例
nu = sum(y_u_r .^ 2) / n2 % Y组原始变量被u_i解释的方差比例
nv = sum(y_v_r .^ 2) / n2 % Y组原始变量被v_i解释的方差比例
fprintf('X组的原始变量被u1-u%d解释的比例为%f\n', num, sum(mu));
fprintf('Y组的原始变量被v1-v%d解释的比例为%f\n', num, sum(nv));

9.2 中国城市竞争力与基础设施的典型相关分析

% 中国城市竞争力与基础设施的典型相关分析
clc; clear
load x.txt
load y.txt
p = size(x, 2); q = size(y, 2);
x = zscore(x); y = zscore(y); % 标准化数据
n = size(x, 1); % 观测数据个数
% 下面做典型相关分析,a1,b1返回的是典型变量的系数,r返回的是典型相关系数
% u1,v1返回的是典型变量的值,stats返回的是假设检验的一些统计量的值
[a1, b1, r, u1, v1, stats] = canoncorr(x, y);
% 下面修正a1,b1的每一列的正负号,使得a,b每一列的系数和为正
% 相应的,典型变量取值的正负号也要修正
a = a1 .* repmat(sign(sum(a1)), size(a1, 1), 1);
b = b1 .* repmat(sign(sum(b1)), size(b1, 1), 1);
u = u1 .* repmat(sign(sum(a1)), size(u1, 1), 1);
v = v1 .* repmat(sign(sum(b1)), size(v1, 1), 1);
x_u_r = x' * u / (n - 1) % 计算x,u的相关系数
y_v_r = y' * v / (n - 1) % 计算y,v的相关系数
x_v_r = x' * v / (n - 1) % 计算x,v的相关系数
y_u_r = y' * u / (n - 1) % 计算y,u的相关系数
ux = sum(x_u_r .^ 2) / p; % X组原始变量被u_i解释的方差比例 
ux_cum = cumsum(ux) % X组原始变量被u_i解释的方差累计比例
vx = sum(x_v_r .^ 2) / p; % X组原始变量被v_i解释的方差比例 
vx_cum = cumsum(vx) % X组原始变量被v_i解释的方差累计比例
vy = sum(y_v_r .^ 2) / q; % Y组原始变量被v_i解释的方差比例 
vy_cum = cumsum(vy) % Y组原始变量被v_i解释的方差累计比例
uy = sum(y_u_r .^ 2) / q; % Y组原始变量被u_i解释的方差比例 
uy_cum = cumsum(uy) % Y组原始变量被u_i解释的方差累计比例
val = r .^ 2 % 典型相关系数的平方,M1或M2矩阵的非零特征值

十、插值与拟合

10.1 机床加工

%% 一维插值
% 一维插值函数 调用格式 y = interp1(x0, y0, x, 'method')
% 其中x0要求单调,method指定插值方法,默认是线性插值,其值为
% nearest:最近项插值 linear:线性插值 spline:三次样条插值 cubic:立方插值
% 机床加工
% x0、y0:题目给出的数据
x0 = [0, 3, 5, 7, 9, 11, 12, 13, 14, 15];
y0 = [0, 1.2, 1.7, 2.0, 2.1, 2.0, 1.8, 1.2, 1.0, 1.6];
x = 0 : 0.1 : 15; % 待插入的数据
y1 = interp1(x0, y0, x); % 待插入数据的纵坐标值
y2 = interp1(x0, y0, x, 'spline');
pp1 = csape(x0, y0);
y3 = fnval(pp1, x);
pp2 = csape(x0, y0, 'second');
y4 = fnval(pp2, x);
subplot(1, 3, 1)
plot(x0, y0, '+', x, y1)
title('Piecewise linear')
subplot(1, 3, 2)
plot(x0, y0, '+', x, y2)
title('Spline1')
subplot(1, 3, 3)
plot(x0, y0, '+', x, y3)
title('Spline2')
dy = diff(y3); dx = diff(x); dy_dx = dy ./ dx;
k0 = dy_dx(1); % x = 0处曲线的斜率
x4 = x(x <= 15 & x >= 13);
y4 = y3(x <= 15 & x >= 13);
[Ymin, Yindex] = min(y4); % 13-15范围内y的最小值
disp([x(Yindex), Ymin, k0])

10.2 求积分

% 求积分
clc; clear
x0 = 0.15 : 0.01 : 0.18;
y0 = [3.5, 1.5, 2.5, 2.8];
pp = csape(x0, y0); % 默认的边界条件是Lagrange边界条件
format long g
xishu = pp.coefs % 显示每个区间上三次多项式的系数
s = integral(@(t)ppval(pp, t), 0.15, 0.18); % 求积分
format % 恢复短小数的显示格式

10.3 丘陵地带

%% 二维插值
% 调用格式：z = interp2(x0, y0, z0, x, y, 'method')
% x0, y0为m维和n维向量,z0是n×m矩阵,表示节点值,x,y为一维数组,表示插值点
% x是行向量,y是列向量,z是矩阵.它的行数是y的维数,列数是x的维数
% method的取值和一维插值法一样。
% 三次样条插值pp = csape({x0, y0}, z0, conds, valconds); z = fnval(pp, {x, y})
% 当插值节点混乱时,使用Zi = griddata(x, y, z, Xi, Yi)
% 丘陵地带
clc; clear
x0 = 100 : 100 : 500;
y0 = 100 : 100 : 400;
z0 = [636, 697, 624, 478, 450;
      698, 712, 630, 478, 420;
      680, 674, 598, 412, 400;
      662, 626, 552, 334, 310];
x = 100 : 10 : 500; y = 100 : 10 : 400;
pp = csape({x0, y0}, z0');
z = fnval(pp, {x, y});
[i, j] = find(z == max(max(z))); % 找最高点的地址 
x1 = x(i); y1 = y(j); zmax = z(i, j); % 求最高点的坐标
% 画三维图
[X, Y] = meshgrid(x, y);
Z = z'; mesh(X, Y, Z); % 不清楚具体用法

10.4 插值节点混乱

%插值节点混乱
x = [129, 140, 103.5, 88, 185.5, 195, 105, 157.5, 107.5, 77,81, 162, 162, 117.5];
y = [7.5, 141.5, 23, 147, 22.5, 137.5, 85.5, -6.5, -81,3, 56.5, -66.5, 84, -33.5];
z = -[4, 8, 6, 8, 6, 8, 8, 9, 9, 8, 8, 9, 4, 9]; % 海底故取负值
% minmax:用于获取数组中每一行最小值和最大值
xmm = minmax(x); ymm = minmax(y);  
X = xmm(1) : xmm(2); Y = ymm(1) : ymm(2);
Z1 = griddata(x, y, z, X, Y', 'cubic'); % 立方插值
Z2 = griddata(x, y, z, X, Y', 'nearest'); % 最近点插值
Z = Z1; % 立方插值和最近点插值的混合插值的初始值
Z(isnan(Z1)) = Z2(isnan(Z1)); % 把立方插值中不确定值换成最近点插值的结果
subplot(1, 2, 1); plot(x, y, '*')
subplot(1, 2, 2); mesh(X, Y, Z)

10.5 曲线拟合

%% 曲线拟合
% 调用格式:a = polyfit(x0, y0, m),x0,y0是数据向量,m表示多项式的阶数
% 返回的是多项式的降幂系数a = [a1, a2, …, am, am+1]
% 要计算x处对应的多项式取值,需调用y = polyval(a, x)
%% 最小二乘曲线拟合
% 调用格式:x = lsqcurvefit(fun, x0, xdata, ydata, lb, ub, options)
% fun:为待拟合的函数,一般定义为F(x,xdata)的M文件
% x0:参数初值值,一般用行向量给出
% xdata,ydata:为观测值,列向量表示
% lb:为参数的下界
% ub:为参数的上界
% options:为选项 
% 解方程组
x = [19, 25, 31, 38, 44]';
y = [19.0, 32.3, 49.0, 73.3, 97.8]';
r = [ones(5, 1), x .^ 2];
ab = r \ y;
x0 = 19 : 0.01 : 44;
y0 = ab(1) + ab(2) * x0 .^ 2;
plot(x, y, 'o', x0, y0, 'r')
% 多项式拟合
x0 = [1990, 1991, 1992, 1993, 1994, 1995, 1996];
y0 = [70, 122, 144, 152, 174, 196, 202];
plot(x0, y0, '*')
a = polyfit(x0, y0, 1);
y97 = polyval(a, 1997);
y98 = polyval(a, 1998);

10.6 黄河小浪底调水调沙问题

% 黄河小浪底调水调沙问题
clc; clear
load data3.txt
liu = data3(1, :)'; % 提出水流量
sha = data3(2, :)'; % 提出含沙量
y = sha .* liu; % 计算排沙量
i = 1 : 24;
t = (12 * i - 4) * 3600;
t1 = t(1); t2 = t(end);
pp = csape(t, y); % 进行三次样条插值
xsh = pp.coefs % 求得插值多项式的系数矩,每一行是一个区间上多项式的系数
TL = integral(@(tt)fnval(pp, tt), t1, t2) % 求总含沙量的积分运算
% 绘制散点图
subplot(1, 2, 1); plot(liu(1 : 11), y(1 : 11), '*') % 第一阶段
subplot(1, 2, 2); plot(liu(12 : 24), y(12 : 24), '*') % 第二阶段
% 进行拟合
format long e
% 以下是第一阶段的拟合
for j = 1 : 2
    nihei1{j} = polyfit(liu(1 : 11), y(1 : 11), j); % 拟合多项式,系数排列从高次幂到低次幂
    yhat1{j} = polyval(nihei1{j}, liu(1 : 11)); % 求预测值
    % 以下求误差平方和与剩余平方差
    cha1(j) = sum((y(1 : 11) - yhat1{j}) .^ 2);
    rmse1(j) = sqrt(cha1(j) / (10 - j));
end
celldisp(nihei1) % 显示细胞数组的所有元素
rmse1
% 以下是第二阶段的拟合
for j = 1 : 2
    nihei2{j} = polyfit(liu(12 : 24), y(12 : 24), (j + 1)); % 这里使用细胞数组
    yhat2{j} = polyval(nihei2{j}, liu(12 : 24)); % 求预测值
    rmse2(j) = sqrt(sum((y(12 : 24) - yhat2{j}) .^ 2) / (11 - j)); % 求剩余标准差
end
celldisp(nihei2) % 显示细胞数组的所有元素
rmse2
format % 恢复默认的短小数的显示格式