最小二乘法拟合出二阶响应面近似模型

背景：根据样本试验数据拟合出二阶响应面近似模型（正交二次型），并使用决定系数R²和调整的决定系数R²_adj来判断二阶响应面模型的拟合精度。

1、样本数据（来源：硕士论文《航空发动机用W形金属密封环密封性能分析与优化》）

编号 x1	     x2	     x3	     x4    	y1	    y2	    y3
1	0.67	0.70	1.78	0.83	804.8	246.7	37.7
2	0.59	0.73	2.09	0.72	686.3	213.5	40.1
3	0.69	0.78	1.81	0.86	703.6	221.0	39.3
4	0.57	0.70	1.84	0.73	822.3	253.6	37.1
5	0.58	0.80	2.07	0.74	632.4	201.6	40.9
6	0.66	0.77	1.95	0.71	716.9	224.7	39.7
7	0.61	0.69	2.04	0.75	730.1	223.3	39.4
8	0.66	0.83	1.88	0.76	676.6	215.8	40.0
9	0.71	0.66	1.85	0.77	845.0	254.6	37.8
10	0.64	0.72	1.73	0.82	813.2	251.8	37.2
11	0.71	0.71	1.99	0.81	709.0	217.4	40.1
12	0.60	0.77	1.91	0.79	698.2	219.5	39.2
13	0.70	0.65	1.94	0.80	794.9	237.9	38.7
14	0.74	0.74	1.96	0.89	663.2	205.7	40.8
15	0.74	0.67	1.74	0.87	839.3	253.8	37.5
16	0.56	0.68	1.87	0.84	773.6	221.6	37.7
17	0.62	0.84	2.06	0.77	598.0	193.0	41.8
18	0.64	0.72	2.01	0.78	703.4	217.3	39.8
19	0.68	0.79	1.79	0.70	777.5	245.8	38.4
20	0.65	0.84	1.81	0.90	640.6	205.3	40.1
21	0.73	0.75	1.76	0.88	743.3	230.9	38.8
22	0.63	0.82	1.70	0.84	726.0	231.4	38.3
23	0.60	0.83	1.90	0.81	648.0	207.5	40.1
24	0.56	0.76	1.97	0.76	693.6	217.9	39.3
25	0.72	0.81	2.03	0.85	600.7	191.7	42.2

2、将上述数据存放在sample_data.xlsx表格中，再编写matlab代码

% % 读取数据（假设数据已经保存在一个名为'data.txt'的文件中，使用空格或制表符分隔）
%% 近似模型：y = a0+a1*x1+a2*x2+a3*x3+a4*x4+a5*x1*x2+a6*x1*x3+a7*x1*x4+a8*x2*x3+a9*x2*x4+a10*x3*x4+a11*x1^2+a12*x2^2+a13*x3^2+a14*x4^2；
clc;clear;
filename = 'sample_data.xlsx';
data = xlsread(filename);

% 提取自变量和因变量
X = data(:, 2:5); % x1, x2, x3, x4
y = data(:, 6:8); % y1, y2, y3

% 构建扩展自变量矩阵（包括线性项、交互项和平方项）
[n, m] = size(X);
X_extended = ones(n, 1); % 初始化扩展矩阵，第一列为截距项

% 添加线性项
X_extended = [X_extended, X];

% 添加交互项
X_extended = [X_extended, X(:,1).*X(:,2), X(:,1).*X(:,3), X(:,1).*X(:,4), ...
              X(:,2).*X(:,3), X(:,2).*X(:,4), X(:,3).*X(:,4)];

% 添加平方项
X_extended = [X_extended, X(:,1).^2, X(:,2).^2, X(:,3).^2, X(:,4).^2];

% 使用最小二乘法求解回归系数（对每个因变量分别进行）
beta_values = zeros(15, 3); % 15个系数（包括截距项），3个因变量
R2_values = zeros(1, 3);
adjusted_R2_values = zeros(1, 3);

for i = 1:3
    % 拟合模型
    beta = (X_extended' * X_extended) \ (X_extended' * y(:, i));
    beta_values(:, i) = beta;
    
    % 预测值
    y_pred = X_extended * beta;
    
    % 计算SST, SSE, SSR
    SST = sum((y(:, i) - mean(y(:, i))).^2);
    SSE = sum((y(:, i) - y_pred).^2);
    SSR = SST - SSE;
    
    % 计算R^2
    R2_values(i) = SSR / SST;
    
    % 计算调整后的R^2
    p = size(X_extended, 2); % 自变量数量（包括截距项）
    adjusted_R2_values(i) = 1 - (1 - R2_values(i)) * ((n - 1) / (n - p));
end

% 显示结果
disp('回归系数（包括截距项）：');
disp(beta_values);
disp('决定系数 R^2:');
disp(R2_values);
disp('调整的决定系数 Adjusted R^2:');
disp(adjusted_R2_values);

3、运行结果