柚子快報(bào)激活碼778899分享：MATLAB算法與模型學(xué)習(xí)筆記

http://yzkb.51969.com/

文章目錄

一、主成分分析二、整數(shù)規(guī)劃2.1 入門題2.2 典型例題2.2.1 背包問題(0-1規(guī)劃)2.2.2 指派問題(0-1規(guī)劃)2.2.3 鋼管切割問題(混合使用切割方案)

三、線性規(guī)劃3.1 入門題3.2 典型例題3.2.1 生產(chǎn)決策問題3.2.2 生成決策問題（整數(shù)規(guī)劃）

四、圖論4.1 最短路徑4.2 最小生成樹4.2.1 prim算法4.2.2 Kruskal算法

4.3 Matlab的圖論工具箱4.4 渡河問題4.5 有向圖最短路徑及長(zhǎng)度4.6 最小生成樹（工具箱）4.7 最大流4.8 旅行商問題

五、聚類分析5.1 DBSCAN算法5.2 聚類分析程序5.3 聚類分析（matlab統(tǒng)計(jì)工具箱）5.4 變量聚類5.5 我國(guó)各地區(qū)普通高等教育發(fā)展情況分析5.5.1 R型聚類(選取指標(biāo))5.5.2 Q型聚類

六、高斯擬合七、非線性規(guī)劃7.1 入門題7.2 典型題7.2.1 選址問題7.2.2 飛行管理問題

八、多目標(biāo)規(guī)劃九、典型相關(guān)分析9.1 職業(yè)滿意度典型相關(guān)分析9.2 中國(guó)城市競(jìng)爭(zhēng)力與基礎(chǔ)設(shè)施的典型相關(guān)分析

十、插值與擬合10.1 機(jī)床加工10.2 求積分10.3 丘陵地帶10.4 插值節(jié)點(diǎn)混亂10.5 曲線擬合10.6 黃河小浪底調(diào)水調(diào)沙問題

一、主成分分析

%% 主成分分析程序

% 主成分回歸分析

clc, clear

load sn.txt

[m, n] = size(sn);

x0 = sn(:, 1 : n-1); y0 = sn(:, n);

hg1 = [ones(m, 1), x0] \ y0; % 計(jì)算普通最小二乘法回歸系數(shù)

hg1 = hg1'; % 變成行向量顯示回歸系數(shù),其中第一個(gè)分量是常數(shù)項(xiàng)

fprintf('y = %f', hg1(1)); % 顯示普通最小二乘法回歸結(jié)果

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); % 計(jì)算相關(guān)系數(shù)矩陣

xd = zscore(x0); % 對(duì)設(shè)計(jì)矩陣進(jìn)行標(biāo)準(zhǔn)化處理

yd = zscore(y0); % 對(duì)y0進(jìn)行標(biāo)準(zhǔn)化處理

[vec1, lamda, rate] = pcacov(r); % vec1特征向量,lamda特征值,rate各主成分的貢獻(xiàn)率

f = repmat(sign(sum(vec1)), size(vec1, 1), 1); % 構(gòu)造與vec1同維度的元素為±1的矩陣

vec2 = vec1 .* f; % 修改特征向量的正負(fù)號(hào),使得特征向量的所有和為正

contr = cumsum(rate); % 計(jì)算累計(jì)貢獻(xiàn)率

df = xd * vec2; % 計(jì)算所有主成分的得分

num = input('請(qǐng)選擇主成分的個(gè)數(shù)：'); % 通過累計(jì)貢獻(xiàn)率交互式選擇主成分的個(gè)數(shù)

hg21 = df(:, 1 : num) \ yd; % 主成分變量的回歸系數(shù),這里由于數(shù)據(jù)標(biāo)準(zhǔn)化,回歸方程的常數(shù)項(xiàng)為0

hg22 = vec2(:, 1 : num) * hg21; % 標(biāo)準(zhǔn)化回歸方程的系數(shù)

% 計(jì)算原始變量回歸方程的系數(shù)

hg23 = [mean(y0) - std(y0) * mean(x0) ./ std(x0) * hg22, std(y0) * hg22' ./ std(x0)];

fprintf('y = %f', hg23(1)); % 顯示主成分回歸結(jié)果

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')

% 下面計(jì)算兩種回歸分析的剩余標(biāo)準(zhǔn)差

rmse1 = sqrt(sum((hg1(1) + x0 * hg1(2 : end)' - y0).^ 2) / (m -n)); % 擬合n個(gè)參數(shù)

rmse2 = sqrt(sum((hg23(1) + x0 * hg23(2 : end)' - y0).^ 2) / (m -num)); % 擬合num個(gè)參數(shù)

二、整數(shù)規(guī)劃

2.1 入門題

clear; clc

%% 第一題

C = [-20, -10]';

intcon = [1, 2]; % 限制第一個(gè)變量和第二個(gè)變量為整數(shù)

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; % 限制第三個(gè)變量為整數(shù)

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規(guī)劃)

C = [-3, -2, -1]';

intcon = 3; % 表示第三個(gè)變量為整數(shù)

A = [1, 1, 1];

b = 7;

Aeq = [4, 2, 1];

beq = 12;

lb = [0, 0, 0]';

ub = [+inf, +inf, 1]'; % 表示0-1規(guī)劃

[x3, fval3] = intlinprog(C, intcon, A, b, Aeq, beq, lb, ub);

2.2 典型例題

2.2.1 背包問題(0-1規(guī)劃)

clear; clc

%% 背包問題(0-1規(guī)劃)

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規(guī)劃)

%% 指派問題(0-1規(guī)劃)

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)]; % 簡(jiǎn)潔

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);

三、線性規(guī)劃

3.1 入門題

clear; clc

%% 第一題

C = [-5, -4, -6]'; % 與標(biāo)準(zhǔn)型保持一致,不轉(zhuǎn)置也行

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]'; % 與標(biāo)準(zhǔn)型保持一致,不轉(zhuǎn)置也行

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]'; % 與標(biāo)準(zhǔn)型保持一致,不轉(zhuǎn)置也行

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; % 結(jié)果取最大值,需添加符號(hào)

%% 多個(gè)解的情況

% 例如:min z = x1 + x2 s.t. x1 + x2 >= 10

C = [1, 1]';

A = [-1, -1];

b = -10;

[x4, fval4] = linprog(C, A, b); % 有多個(gè)解時(shí),返回一個(gè)值

%% 不存在解的情況

% 例如: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 生產(chǎn)決策問題

clear; clc

%% 生產(chǎn)決策問題(本題應(yīng)為整數(shù)規(guī)劃,僅供入門參考)

format long g % 計(jì)算結(jié)果顯示為長(zhǎng)數(shù)字

% 目標(biāo)函數(shù)系數(shù)向量

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;

% 系數(shù)矩陣(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;

% 常數(shù)項(xiàng)(b) 不等式約束

b = zeros(5, 1);

b(1) = 1200;

b(2) = 10000;

b(3) = 2000;

b(4) = 7000;

b(5) = 4000;

% 系數(shù)矩陣(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;

% 常數(shù)項(xiàng)(beq) 等式約束

beq = zeros(2, 1);

% 下界(lb)

lb = zeros(9, 1);

[x1, fval1] = linprog(-C, A, b, Aeq, beq, lb);

fval1 = -fval1; % 題目為求最大值

3.2.2 生成決策問題（整數(shù)規(guī)劃）

% 將其改為整數(shù)規(guī)劃

intcon = 1 : 9;

[x, fval] = intlinprog(-C, intcon, A, b, Aeq, beq, lb);

fval = -fval;

%% 投料問題(運(yùn)輸問題)

format long g % 計(jì)算結(jié)果顯示為長(zhǎng)數(shù)字

% 目標(biāo)函數(shù)系數(shù)向量

M = [5, 1;

2, 7]; % 產(chǎn)地

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

% 系數(shù)矩陣(A) 不等式約束

A = zeros(2, 12);

A(1, 1 : 6) = 1; % 簡(jiǎn)潔

A(2, 7 : 12) = 1; % 簡(jiǎn)潔

% 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;

% 常數(shù)項(xiàng)(b) 不等式約束

b = [20, 20]';

% 系數(shù)矩陣(Aeq) 等式約束

Aeq = zeros(6, 12);

for i = 1 : 6

Aeq(i, i) = 1; Aeq(i, i + 6) = 1;

end

Aeq = [eye(6), eye(6)]; % 簡(jiǎn)潔

% 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;

% 常數(shù)項(xiàng)(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算法:求圖中任意兩點(diǎn)之間最短路徑

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); % 找第一個(gè)最小的值

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的圖論工具箱

% 無向圖最短路徑及長(zhǎng)度

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工具箱要求數(shù)據(jù)是下三角矩陣

[i, j, v] = find(a);

b = sparse(i, j, v, 11, 11); % 構(gòu)造系數(shù)矩陣

% Directed是標(biāo)志圖為有向或無向的屬性

[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]; % 每一行是一個(gè)可行狀態(tài)

b = [1, 0, 0, 0;

1, 1, 0, 0;

1, 0, 1, 0;

1, 0, 0, 1]; % 每一行是一個(gè)可以轉(zhuǎn)移狀態(tài)

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); % 構(gòu)造稀疏矩陣

% 該圖是無向圖,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); % 線型寬度設(shè)置為1.5

4.5 有向圖最短路徑及長(zhǎng)度

% 有向圖最短路徑及長(zhǎng)度

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); % 構(gòu)造稀疏矩陣

% 有向圖Directed屬性值為真或1,方法屬性的默認(rèn)值是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的兩兩列向量間的絕對(duì)值距離

d = tril(d); % 截取Matlab工具箱要求的下三角形矩陣

b = sparse(d); % 轉(zhuǎn)化為稀疏矩陣

[ST, pred] = graphminspantree(b, 'Method', 'Kruskal') % 調(diào)用最小生成樹的命令

st = full(ST); % 把最小生成樹的稀疏矩陣轉(zhuǎn)化為普通矩陣

TreeLength = sum(sum(st)) % 求最小生成樹長(zhǎng)度

view(biograph(ST, [], 'ShowArrows', 'off')) % 畫出最小生成樹

4.7 最大流

% 最大流

% 只能解決權(quán)重都為正值,且兩個(gè)頂點(diǎn)之間不能有兩條弧

% 添加虛擬頂點(diǎn)9,解決頂點(diǎn)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; % 退出標(biāo)志

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; % 修改一次,標(biāo)志加1

end

end

end

if flag == 0 % 一條邊也沒有修改,就返回

long = 0; % 圈長(zhǎng)的初始值

for i = 1 : L

long = long + a(c(i), c(i + 1)); % 求改良圈的長(zhǎng)度

end

circle = c; % 返回修改圈

return

end

end

五、聚類分析

5.1 DBSCAN算法

%% DBSCAN算法

clc; clear;

load mydata; % 加載數(shù)據(jù)

epsilon = 0.5; % 半徑

MinPts = 10; % 點(diǎn)數(shù)

C = 0; % 記錄分的類數(shù)

n = size(X, 1); % X為加載的數(shù)據(jù)集

IDX = zeros(n, 1); % 初始化全部為0,即全部為噪音點(diǎn)

D = pdist2(X, X); % 計(jì)算每個(gè)點(diǎn)之間的距離

visited = false(n, 1); % 是否已經(jīng)訪問

isnoise = false(n, 1); % 是否為噪音點(diǎn)

for i = 1 : n

if ~visited(i)

visited(i) = true; % 更新為已訪問

% 記錄該訪問點(diǎn)的臨近點(diǎn)

Neighbors = find(D(i, :) <= epsilon);

if numel(Neighbors) < MinPts % numel()計(jì)算數(shù)組中元素?cái)?shù)量

% 附近點(diǎn)的數(shù)量小于已給數(shù)量,記為噪聲點(diǎn)

isnoise(i) = true;

else

C = C + 1;

IDX(i) = C;

% 開始搜索其臨近點(diǎn)

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

% 如果只有兩個(gè)指標(biāo)的,就可以繪制圖形

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 % 噪聲點(diǎn)

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; % 增加網(wǎng)格線

legend(Legends); % 標(biāo)簽

5.2 聚類分析程序

%% 聚類分析程序

clc, clear

a = [1, 0; 1, 1; 3, 2; 4, 3; 2, 5];

[m, n] = size(a);

d = mandist(a'); % mandist求矩陣列向量組之間的兩兩絕對(duì)值距離

d = tril(d); % 截取下三角形

nd = nonzeros(d); % 去掉d中的零元素,非零元素按列排列

d = union([], nd); % 去掉重復(fù)的非零元素

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)); % 把數(shù)組tm變成行向量

fprintf('第%d次合成,平臺(tái)高度為%d時(shí)的分類結(jié)果為：%s\n', i, nd_min, int2str(tm));

nd(nd == nd_min) = []; % 刪除已經(jīng)歸類的元素

if length(nd) == 1

break

end

end

5.3 聚類分析（matlab統(tǒng)計(jì)工具箱）

% 使用matlab統(tǒng)計(jì)工具箱

clc, clear

a = [1, 0; 1, 1; 3, 2; 4, 3; 2, 5];

y = pdist(a, 'cityblock'); % 求a的兩兩行向量間的絕對(duì)值距離

yc = squareform(y); % 變換成距離矩陣

z = linkage(y); % 產(chǎn)生等級(jí)聚類書

dendrogram(z) % 畫聚類圖

T = cluster(z, 'maxclust', 3); % 把對(duì)象劃分成3類

% zsore(X) 對(duì)數(shù)據(jù)矩陣進(jìn)行標(biāo)準(zhǔn)化處理

for i = 1 : 3

tm = find(T == i); % 求第i類的對(duì)象

tm = reshape(tm, 1, length(tm)); % 變成行向量

fprintf('第%d類的有%s\n', i, int2str(tm)); %顯示分類結(jié)果

end

5.4 變量聚類

% 對(duì)變量聚類

clc, clear

a = textread('ch.txt');

d = 1 - abs(a); % 進(jìn)行數(shù)據(jù)變換,把相關(guān)系數(shù)轉(zhuǎn)化為距離

d = tril(d); % 提取出d矩陣的下三角部分

b = nonzeros(d); % 去掉d中的零元素

b = b'; % 化為行向量(對(duì)變量聚類)

z = linkage(b, 'complete'); % 按最長(zhǎng)距離法聚類

y = cluster(z, 'maxclust', 2); % 把變量劃分為兩類

ind1 = find(y == 1); ind1 = ind1'; % 顯示第一類對(duì)應(yīng)的變量標(biāo)號(hào)

ind2 = find(y == 2); ind2 = ind2'; % 顯示第二類對(duì)應(yīng)的變量標(biāo)號(hào)

h = dendrogram(z); % 畫聚類圖

set(h, 'Color', 'k', 'LineWidth', 1.3) % 把聚類圖線的顏色改為黑色,線寬加粗

5.5 我國(guó)各地區(qū)普通高等教育發(fā)展情況分析

5.5.1 R型聚類(選取指標(biāo))

% R型聚類(選取指標(biāo))

clc, clear

a = load('gj.txt'); % 把原始數(shù)據(jù)保存在純文本文件gj.txt中

b = zscore(a); % 數(shù)據(jù)標(biāo)準(zhǔn)化

r = corrcoef(b); % 計(jì)算相關(guān)系數(shù)矩陣

% d = tril(1 - r); d = nonzeros(d)'; % 另一種計(jì)算距離方法

d = pdist(b', 'correlation'); % 計(jì)算相關(guān)系數(shù)到處的距離

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類對(duì)象

tm = reshape(tm, 1, length(tm)); % 變成行向量

fprintf('第%d類的有%s\n', i, int2str(tm)); % 顯示分類結(jié)果

end

5.5.2 Q型聚類

% Q型聚類

clc, clear

load gj.txt % 把原始數(shù)據(jù)保存在純文本文件gj.txt中

gj(:, 3 : 6) = []; % 刪除數(shù)據(jù)矩陣的第3-6列,即使用變量1,2,7,8,9,10

gj = zscore(gj); % 數(shù)據(jù)標(biāo)準(zhǔn)化

y = pdist(gj); % 求對(duì)象間的歐氏距離,每行是一個(gè)對(duì)象

z = linkage(y, 'average'); % 按類平均法聚類

h = dendrogram(z); % 畫聚類圖

set(h, 'Color', 'k', 'LineWidth', 1.3) % 把聚類圖線的顏色改為黑色,線寬加粗

for k = 3 : 5

fprintf('劃分為%d類的結(jié)果如下：\n', k)

T = cluster(z, 'maxclust', k); % 把樣本點(diǎn)劃分為k類

for i = 1 : k

tm = find(T == i); % 求第i類對(duì)象

tm = reshape(tm, 1, length(tm)); % 變成行向量

fprintf('第%d類的有%s\n', i, int2str(tm)); % 顯示分類結(jié)果

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; % 定義區(qū)間的填充顏色

h1.EdgeColor = [1, 1, 1]; % 邊界顏色設(shè)置為白色

alpha .2 % 設(shè)置透明色

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')

七、非線性規(guī)劃

7.1 入門題

%% 非線性規(guī)劃

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算法(內(nèi)點(diǎn)法)

option = optimoptions('fmincon', 'Algorithm', 'interior-point');

[x2, fval2] = fmincon(@fun1, x0, A, b, [], [], [], [], @nonlfun1, option);

fval2 = -fval2;

% 使用SQP算法(序列二次規(guī)劃法)

option = optimoptions('fmincon', 'Algorithm', 'sqp');

[x3, fval3] = fmincon(@fun1, x0, A, b, [], [], [], [], @nonlfun1, option);

fval3 = -fval3; % 結(jié)果較差

x00 = [1, 1]; % 更改初始值

[x4, fval4] = fmincon(@fun1, x00, A, b, [], [], [], [], @nonlfun1, option);

fval4 = -fval4; % 效果變好,說明序列二次規(guī)劃法適應(yīng)性較差

% 使用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

% 系數(shù)矩陣(A) 不等式約束

A = zeros(2, 16);

A(1, 1 : 6) = 1; % 簡(jiǎn)潔

A(2, 7 : 12) = 1; % 簡(jiǎn)潔

% 常數(shù)項(xiàng)(b) 不等式約束

b = [20, 20]';

% 系數(shù)矩陣(Aeq) 等式約束

Aeq = [eye(6), eye(6), zeros(6, 4)]; % 簡(jiǎn)潔

% 常數(shù)項(xiàng)(beq) 等式約束

beq = [3, 5, 4, 7, 6, 11]';

% 下界(lb)

lb = zeros(16, 1);

% 進(jìn)行求解

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

% 畫六架飛機(jī)的位置

figure(1) % 生成一個(gè)圖形

box on % 顯示為封閉的盒子

% 繪制飛機(jī)的初始位置

data = [150, 140, 243;

85, 85, 236;

150, 155, 220.5;

145, 50, 159;

130, 150, 230;

0, 0, 52]; % x坐標(biāo) y坐標(biāo) 方向角

plot(data(:, 1), data(:, 2), '.r')

axis([0, 160, 0, 160]) % 設(shè)置坐標(biāo)軸刻度范圍

hold on % 接著繪制

% 在圖上標(biāo)上注釋

for i = 1 : 6

txt = ['飛機(jī)', num2str(i)];

text(data(i, 1) + 2, data(i, 2) + 2, txt, 'FontSize', 8)

end

% 求解非線性規(guī)劃問題

x0 = [0, 0, 0, 0, 0, 0]; % 飛機(jī)調(diào)整角度的初始值

lb = -pi / 6 * ones(6, 1); % 下界

ub = pi / 6 * ones(6, 1); % 上界

[x2, fval2] = fmincon(@fun5, x0, [], [], [], [], lb, ub, @nonlfun5);

x2 = x2 * pi / 180; % 將弧度制轉(zhuǎn)換為角度制

八、多目標(biāo)規(guī)劃

%% 第一種方式

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]';

% 第一個(gè)目標(biāo)函數(shù)的目標(biāo)值

[x1, g1] = linprog(c1, a, b, [], [], zeros(4, 1));

% 第二個(gè)目標(biāo)函數(shù)的目標(biāo)值

[x2, g2] = linprog(c2, a, b, [], [], zeros(4, 1));

g3 = [g1, g2]'; % 目標(biāo)goal的值

% 這里的權(quán)重weight=目標(biāo)goal的絕對(duì)值

[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; % 用匿名函數(shù)定義目標(biāo)向量

% 第一個(gè)目標(biāo)函數(shù)的目標(biāo)值

[x1, g1] = linprog(c1, a, b, [], [], zeros(4, 1));

% 第二個(gè)目標(biāo)函數(shù)的目標(biāo)值

[x2, g2] = linprog(c2, a, b, [], [], zeros(4, 1));

g3 = [g1, g2]'; % 目標(biāo)goal的值

% 這里的權(quán)重weight=目標(biāo)goal的絕對(duì)值

[x, fval] = fgoalattain(Fun, rand(4, 1), g3, abs(g3), a, b, [], [], zeros(4, 1));

九、典型相關(guān)分析

9.1 職業(yè)滿意度典型相關(guān)分析

%% 典型相關(guān)分析

% 職業(yè)滿意度典型相關(guān)分析

clc; clear

load r.txt

n1 = 5; n2 = 7; num = min(n1, n2);

s11 = r(1 : n1, 1 : n1); % 提出X與X的相關(guān)系數(shù)

s12 = r(1 : n1, n1 + 1 : end); % 提出X與Y的相關(guān)系數(shù)

s21 = s12'; % 提出Y與X的相關(guān)系數(shù)

s22 = r(n1 + 1 : end, n1 + 1 : end); % 提出Y與Y的相關(guān)系數(shù)

m1 = inv(s11) * s12 * inv(s22) * s21; % 計(jì)算矩陣M1

m2 = inv(s22) * s21 * inv(s11) * s12; % 計(jì)算矩陣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)); % 計(jì)算特征值的平方根

[val1, ind1] = sort(val1, 'descend'); % 按照從大到小排列

a = vec1(:, ind1(1 : num)); % 取出X組的系數(shù)陣

dcoef1 = val1(1 : num); % 提出典型相關(guān)系數(shù)

xlswrite('bk.xls', a, 'Sheet1', 'A1') % 把計(jì)算結(jié)果寫到Excel文件中去

flag = n1 + 2; % % 把計(jì)算結(jié)果寫到Excel中的行計(jì)數(shù)變量

str = char(['A', int2str(flag)]); % str為Excel中寫數(shù)據(jù)的起始位置

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)); % 計(jì)算特征值的平方根

[val2, ind2] = sort(val2, 'descend'); % 按照從大到小排列

b = vec2(:, ind2(1 : num)); % 取出Y組的系數(shù)陣

dcoef2 = val2(1 : num); % 提出典型相關(guān)系數(shù)

flag = flag + 2; str = char(['A', int2str(flag)]); % str為Excel中寫數(shù)據(jù)的起始位置

xlswrite('bk.xls', b, 'Sheet1', str) % 把計(jì)算結(jié)果寫到Excel文件中去

flag = flag + n2 + 1; str = char(['A', int2str(flag)]); % str為Excel中寫數(shù)據(jù)的起始位置

xlswrite('bk.xls', dcoef2', 'Sheet1', str)

x_u_r = s11 * a; % x,u的相關(guān)系數(shù)

y_v_r = s22 * b; % y,v的相關(guān)系數(shù)

x_v_r = s12 * b; % x,v的相關(guān)系數(shù)

y_u_r = s21 * a; % y,u的相關(guān)系數(shù)

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 中國(guó)城市競(jìng)爭(zhēng)力與基礎(chǔ)設(shè)施的典型相關(guān)分析

% 中國(guó)城市競(jìng)爭(zhēng)力與基礎(chǔ)設(shè)施的典型相關(guān)分析

clc; clear

load x.txt

load y.txt

p = size(x, 2); q = size(y, 2);

x = zscore(x); y = zscore(y); % 標(biāo)準(zhǔn)化數(shù)據(jù)

n = size(x, 1); % 觀測(cè)數(shù)據(jù)個(gè)數(shù)

% 下面做典型相關(guān)分析,a1,b1返回的是典型變量的系數(shù),r返回的是典型相關(guān)系數(shù)

% u1,v1返回的是典型變量的值,stats返回的是假設(shè)檢驗(yàn)的一些統(tǒng)計(jì)量的值

[a1, b1, r, u1, v1, stats] = canoncorr(x, y);

% 下面修正a1,b1的每一列的正負(fù)號(hào),使得a,b每一列的系數(shù)和為正

% 相應(yīng)的,典型變量取值的正負(fù)號(hào)也要修正

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) % 計(jì)算x,u的相關(guān)系數(shù)

y_v_r = y' * v / (n - 1) % 計(jì)算y,v的相關(guān)系數(shù)

x_v_r = x' * v / (n - 1) % 計(jì)算x,v的相關(guān)系數(shù)

y_u_r = y' * u / (n - 1) % 計(jì)算y,u的相關(guān)系數(shù)

ux = sum(x_u_r .^ 2) / p; % X組原始變量被u_i解釋的方差比例

ux_cum = cumsum(ux) % X組原始變量被u_i解釋的方差累計(jì)比例

vx = sum(x_v_r .^ 2) / p; % X組原始變量被v_i解釋的方差比例

vx_cum = cumsum(vx) % X組原始變量被v_i解釋的方差累計(jì)比例

vy = sum(y_v_r .^ 2) / q; % Y組原始變量被v_i解釋的方差比例

vy_cum = cumsum(vy) % Y組原始變量被v_i解釋的方差累計(jì)比例

uy = sum(y_u_r .^ 2) / q; % Y組原始變量被u_i解釋的方差比例

uy_cum = cumsum(uy) % Y組原始變量被u_i解釋的方差累計(jì)比例

val = r .^ 2 % 典型相關(guān)系數(shù)的平方,M1或M2矩陣的非零特征值

十、插值與擬合

10.1 機(jī)床加工

%% 一維插值

% 一維插值函數(shù) 調(diào)用格式 y = interp1(x0, y0, x, 'method')

% 其中x0要求單調(diào),method指定插值方法,默認(rèn)是線性插值,其值為

% nearest:最近項(xiàng)插值 linear:線性插值 spline:三次樣條插值 cubic:立方插值

% 機(jī)床加工

% x0、y0:題目給出的數(shù)據(jù)

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; % 待插入的數(shù)據(jù)

y1 = interp1(x0, y0, x); % 待插入數(shù)據(jù)的縱坐標(biāo)值

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范圍內(nèi)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); % 默認(rèn)的邊界條件是Lagrange邊界條件

format long g

xishu = pp.coefs % 顯示每個(gè)區(qū)間上三次多項(xiàng)式的系數(shù)

s = integral(@(t)ppval(pp, t), 0.15, 0.18); % 求積分

format % 恢復(fù)短小數(shù)的顯示格式

10.3 丘陵地帶

%% 二維插值

% 調(diào)用格式：z = interp2(x0, y0, z0, x, y, 'method')

% x0, y0為m維和n維向量,z0是n×m矩陣,表示節(jié)點(diǎn)值,x,y為一維數(shù)組,表示插值點(diǎn)

% x是行向量,y是列向量,z是矩陣.它的行數(shù)是y的維數(shù),列數(shù)是x的維數(shù)

% method的取值和一維插值法一樣。

% 三次樣條插值pp = csape({x0, y0}, z0, conds, valconds); z = fnval(pp, {x, y})

% 當(dāng)插值節(jié)點(diǎn)混亂時(shí),使用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))); % 找最高點(diǎn)的地址

x1 = x(i); y1 = y(j); zmax = z(i, j); % 求最高點(diǎn)的坐標(biāo)

% 畫三維圖

[X, Y] = meshgrid(x, y);

Z = z'; mesh(X, Y, Z); % 不清楚具體用法

10.4 插值節(jié)點(diǎn)混亂

%插值節(jié)點(diǎn)混亂

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]; % 海底故取負(fù)值

% minmax:用于獲取數(shù)組中每一行最小值和最大值

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'); % 最近點(diǎn)插值

Z = Z1; % 立方插值和最近點(diǎn)插值的混合插值的初始值

Z(isnan(Z1)) = Z2(isnan(Z1)); % 把立方插值中不確定值換成最近點(diǎn)插值的結(jié)果

subplot(1, 2, 1); plot(x, y, '*')

subplot(1, 2, 2); mesh(X, Y, Z)

10.5 曲線擬合

%% 曲線擬合

% 調(diào)用格式:a = polyfit(x0, y0, m),x0,y0是數(shù)據(jù)向量,m表示多項(xiàng)式的階數(shù)

% 返回的是多項(xiàng)式的降冪系數(shù)a = [a1, a2, …, am, am+1]

% 要計(jì)算x處對(duì)應(yīng)的多項(xiàng)式取值,需調(diào)用y = polyval(a, x)

%% 最小二乘曲線擬合

% 調(diào)用格式:x = lsqcurvefit(fun, x0, xdata, ydata, lb, ub, options)

% fun:為待擬合的函數(shù),一般定義為F(x,xdata)的M文件

% x0:參數(shù)初值值,一般用行向量給出

% xdata,ydata:為觀測(cè)值,列向量表示

% lb:為參數(shù)的下界

% ub:為參數(shù)的上界

% options:為選項(xiàng)

% 解方程組

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')

% 多項(xiàng)式擬合

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 黃河小浪底調(diào)水調(diào)沙問題

% 黃河小浪底調(diào)水調(diào)沙問題

clc; clear

load data3.txt

liu = data3(1, :)'; % 提出水流量

sha = data3(2, :)'; % 提出含沙量

y = sha .* liu; % 計(jì)算排沙量

i = 1 : 24;

t = (12 * i - 4) * 3600;

t1 = t(1); t2 = t(end);

pp = csape(t, y); % 進(jìn)行三次樣條插值

xsh = pp.coefs % 求得插值多項(xiàng)式的系數(shù)矩,每一行是一個(gè)區(qū)間上多項(xiàng)式的系數(shù)

TL = integral(@(tt)fnval(pp, tt), t1, t2) % 求總含沙量的積分運(yùn)算

% 繪制散點(diǎn)圖

subplot(1, 2, 1); plot(liu(1 : 11), y(1 : 11), '*') % 第一階段

subplot(1, 2, 2); plot(liu(12 : 24), y(12 : 24), '*') % 第二階段

% 進(jìn)行擬合

format long e

% 以下是第一階段的擬合

for j = 1 : 2

nihei1{j} = polyfit(liu(1 : 11), y(1 : 11), j); % 擬合多項(xiàng)式,系數(shù)排列從高次冪到低次冪

yhat1{j} = polyval(nihei1{j}, liu(1 : 11)); % 求預(yù)測(cè)值

% 以下求誤差平方和與剩余平方差

cha1(j) = sum((y(1 : 11) - yhat1{j}) .^ 2);

rmse1(j) = sqrt(cha1(j) / (10 - j));

end

celldisp(nihei1) % 顯示細(xì)胞數(shù)組的所有元素

rmse1

% 以下是第二階段的擬合

for j = 1 : 2

nihei2{j} = polyfit(liu(12 : 24), y(12 : 24), (j + 1)); % 這里使用細(xì)胞數(shù)組

yhat2{j} = polyval(nihei2{j}, liu(12 : 24)); % 求預(yù)測(cè)值

rmse2(j) = sqrt(sum((y(12 : 24) - yhat2{j}) .^ 2) / (11 - j)); % 求剩余標(biāo)準(zhǔn)差

end

celldisp(nihei2) % 顯示細(xì)胞數(shù)組的所有元素

rmse2

format % 恢復(fù)默認(rèn)的短小數(shù)的顯示格式

柚子快報(bào)激活碼778899分享：MATLAB算法與模型學(xué)習(xí)筆記

http://yzkb.51969.com/

文章鏈接