柚子快報激活碼778899分享：MATLAB算法與模型學習筆記

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 有向圖最短路徑及長度4.6 最小生成樹（工具箱）4.7 最大流4.8 旅行商問題

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

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

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

十、插值與擬合10.1 機床加工10.2 求積分10.3 丘陵地帶10.4 插值節(jié)點混亂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; % 計算普通最小二乘法回歸系數(shù)

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

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

xd = zscore(x0); % 對設(shè)計矩陣進行標準化處理

yd = zscore(y0); % 對y0進行標準化處理

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

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

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

contr = cumsum(rate); % 計算累計貢獻率

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

num = input('請選擇主成分的個數(shù)：'); % 通過累計貢獻率交互式選擇主成分的個數(shù)

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

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

% 計算原始變量回歸方程的系數(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')

% 下面計算兩種回歸分析的剩余標準差

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

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

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

2.1 入門題

clear; clc

%% 第一題

C = [-20, -10]';

intcon = [1, 2]; % 限制第一個變量和第二個變量為整數(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; % 限制第三個變量為整數(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; % 表示第三個變量為整數(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)]; % 簡潔

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]'; % 與標準型保持一致,不轉(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]'; % 與標準型保持一致,不轉(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]'; % 與標準型保持一致,不轉(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é)果取最大值,需添加符號

%% 多個解的情況

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

clear; clc

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

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

% 目標函數(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ù)項(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ù)項(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;

%% 投料問題(運輸問題)

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

% 目標函數(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; % 簡潔

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;

% 常數(shù)項(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)]; % 簡潔

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

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

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

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

b = [1, 0, 0, 0;

1, 1, 0, 0;

1, 0, 1, 0;

1, 0, 0, 1]; % 每一行是一個可以轉(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 有向圖最短路徑及長度

% 有向圖最短路徑及長度

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,方法屬性的默認值是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); % 轉(zhuǎn)化為稀疏矩陣

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

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

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

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

4.7 最大流

% 最大流

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

% 添加虛擬頂點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; % 加載數(shù)據(jù)

epsilon = 0.5; % 半徑

MinPts = 10; % 點數(shù)

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

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

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

D = pdist2(X, X); % 計算每個點之間的距離

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

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

for i = 1 : n

if ~visited(i)

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

% 記錄該訪問點的臨近點

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

if numel(Neighbors) < MinPts % numel()計算數(shù)組中元素數(shù)量

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

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

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

fprintf('第%d次合成,平臺高度為%d時的分類結(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)計工具箱）

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

clc, clear

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

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

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

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

dendrogram(z) % 畫聚類圖

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

% zsore(X) 對數(shù)據(jù)矩陣進行標準化處理

for i = 1 : 3

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

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

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

end

5.4 變量聚類

% 對變量聚類

clc, clear

a = textread('ch.txt');

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

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'; % 顯示第一類對應(yīng)的變量標號

ind2 = find(y == 2); ind2 = ind2'; % 顯示第二類對應(yīng)的變量標號

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

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

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

5.5.1 R型聚類(選取指標)

% R型聚類(選取指標)

clc, clear

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

b = zscore(a); % 數(shù)據(jù)標準化

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

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

d = pdist(b', 'correlation'); % 計算相關(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類對象

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ù)標準化

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

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); % 把樣本點劃分為k類

for i = 1 : k

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

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)點法)

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

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

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

b = [20, 20]';

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

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

% 常數(shù)項(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]) % 設(shè)置坐標軸刻度范圍

hold on % 接著繪制

% 在圖上標上注釋

for i = 1 : 6

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

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

end

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

x0 = [0, 0, 0, 0, 0, 0]; % 飛機調(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)換為角度制

八、多目標規(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]';

% 第一個目標函數(shù)的目標值

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

% 第二個目標函數(shù)的目標值

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

g3 = [g1, g2]'; % 目標goal的值

% 這里的權(quán)重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; % 用匿名函數(shù)定義目標向量

% 第一個目標函數(shù)的目標值

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

% 第二個目標函數(shù)的目標值

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

g3 = [g1, g2]'; % 目標goal的值

% 這里的權(quán)重weight=目標goal的絕對值

[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; % 計算矩陣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組的系數(shù)陣

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

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

flag = n1 + 2; % % 把計算結(jié)果寫到Excel中的行計數(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)); % 計算特征值的平方根

[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) % 把計算結(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 中國城市競爭力與基礎(chǔ)設(shè)施的典型相關(guān)分析

% 中國城市競爭力與基礎(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); % 標準化數(shù)據(jù)

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

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

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

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

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

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

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的相關(guān)系數(shù)

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

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

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

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 % 典型相關(guān)系數(shù)的平方,M1或M2矩陣的非零特征值

十、插值與擬合

10.1 機床加工

%% 一維插值

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

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

% nearest:最近項插值 linear:線性插值 spline:三次樣條插值 cubic:立方插值

% 機床加工

% 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ù)的縱坐標值

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); % 默認的邊界條件是Lagrange邊界條件

format long g

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

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

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

10.3 丘陵地帶

%% 二維插值

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

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

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

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

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

% 當插值節(jié)點混亂時,使用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 插值節(jié)點混亂

%插值節(jié)點混亂

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:用于獲取數(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'); % 最近點插值

Z = Z1; % 立方插值和最近點插值的混合插值的初始值

Z(isnan(Z1)) = Z2(isnan(Z1)); % 把立方插值中不確定值換成最近點插值的結(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表示多項式的階數(shù)

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

% 要計算x處對應(yī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:為觀測值,列向量表示

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

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

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

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

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 % 求得插值多項式的系數(shù)矩,每一行是一個區(qū)間上多項式的系數(shù)

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); % 擬合多項式,系數(shù)排列從高次冪到低次冪

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

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

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

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

end

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

rmse1

% 以下是第二階段的擬合

for j = 1 : 2

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

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

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

end

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

rmse2

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

柚子快報激活碼778899分享：MATLAB算法與模型學習筆記

http://yzkb.51969.com/

文章鏈接