clear;clc;

pz = load('pz.txt');

mu = mean(pz);       % 均值

sig = std(pz);           % 方差

rr = corrcoef(pz);      % 相关系数

data = zscore(pz);     % 数据标准化

% data = (pz - mu(ones(20, 1), :))./sig(ones(20, 1), :);

n = 3;                        % 自变量个数

m = 3;                       % 因变量个数

x0 = pz(:, 1:n);

y0 = pz(:, n + 1:end);

e0 = data(:, 1:n);

f0 = data(:, n + 1:end);

num = size(e0, 1);      % 样本点个数

chg = eye(n);          % w到w*变换矩阵的初始化

for i = 1:n

%%%  以下计算w，w*和t的得分向量

matrix = e0'*(f0*f0')*e0;

[vec, val] = eig(matrix);     % 特征值和特征向量

val = diag(val);

[val, ind] = sort(val, 'descend');

w(:, i) = vec(:, ind(1));   % 提取最大特征值和对应的特征向量

w_star(:, i) = chg*w(:, i);                   % 计算w*的取值

t(:, i) = e0*w(:, i);                               % 计算成分ti的得分

alpha = e0'*t(:, i)/(t(:, i)'*t(:, i));           % 计算alpha_i

chg = chg*(eye(n) - w(:, i)*alpha');      % 计算w到w*的变换矩阵

e = e0 - t(:, i)*alpha';                            % 计算残差

e0 = e;

%%% 以下计算ss(i)的值

beta = [t(:, 1:i), ones(num, 1)]\f0;            % 求回归方程的系数

beta(end, :) = [];                                      % 删除回归分析的常数项

cancha = f0 - t(:, 1:i)*beta;                     % 残差

ss(i) = norm(cancha)^2;                           % 误差平方和

%%% 以下计算press(i)

for j = 1:num

t1 = t(:, 1:i);

f1 = f0;

she_t = t1(j, :);              % 保存舍去的第j个样本点

she_f = f1(j, :);

t1(j, :) = [];                     % 删除第j个观测值

f1(j, :) = [];

beta1 = [t1, ones(num - 1, 1)]\f1;        % 求回归分析的系数

beta1(end, :) = [];              % 删除回归分析的常数项

cancha = she_f - she_t*beta1;          % 残差

press_i(j) = norm(cancha)^2;

end

press(i) = sum(press_i);

if i > 1

Q_h2(i) = 1 - press(i)/ss(i - 1);

else

Q_h2(i) = 1;

end

if Q_h2(i) < 0.0975

fprintf('提出的成分个数r=%d', i);

r = i;

break;

end

end

beta_z = [t(:, 1:r), ones(num, 1)]\f0;            % 求Y关于t的回归系数

beta_z(end, :) = [];                                       % 删除常数项

xishu = w_star(:, 1:r)*beta_z;     % 求Y关于标准数据X的回归系数，每一列是一个回归方程

mu_x = mu(1:n);

mu_y = mu(n + 1:end);

sig_x = sig(1:n);

sig_y = sig(n + 1:end);

for i = 1:m

ch0(i) = mu_y(i) - mu_x./sig_x*sig_y(i)*xishu(:, i);     % 计算原始数据的回归方程的常数项

end

for i = 1:m

xish(:, i) = xishu(:, i)./sig_x'*sig_y(i);    % 计算原始数据的回归方程的系数，每一列是一个回归方程

end

sol = [ch0; xish];    % 回归方程的系数，每一列是一个方程，每一列的第一个数是常数项

disp('回归方程的系数，每一列是一个方程，每一列的第一个数是常数项');

disp(sol);

save mydata x0 y0 num xishu ch0 xish

%%%

subplot(2, 2, 1);

bar(xishu');

ch0 = repmat(ch0, num, 1);

yhat = ch0 + x0*xish; %计算y 的预测值

y1max = max(yhat);

y2max = max(y0);

ymax = max([y1max; y2max]);

cancha = yhat - y0; %计算残差

subplot(2, 2, 2)

plot(0:ymax(1), 0:ymax(1), yhat(:, 1), y0(:, 1), '*')

subplot(2, 2, 3)

plot(0:ymax(2), 0:ymax(2), yhat(:,2), y0(:,2), 'O')

subplot(2, 2, 4)

plot(0:ymax(3), 0:ymax(3), yhat(:,3), y0(:,3), 'H')

pz.txt内容

191 36 50   5 162   60

189 37 52   2 110   60

193 38 58 12 101 101

162 35 62 12 105   37

189 35 46 13 155   58

182 36 56   4 101   42

211 38 56   8 101   38

167 34 60   6 125   40

176 31 74 15 200   40

154 33 56 17 251 250

169 34 50 17 120   38

166 33 52 13 210 115

154 34 64 14 215 105

247 46 50   1   50   50

193 36 46   6   70   31

202 37 62 12 210 120

176 37 54   4   60   25

157 32 52 11 230   80

156 33 54 15 225   73

138 33 68   2 110   43

结果：

untitled.png

(18.68 KB, 下载次数: 9)

2017-5-16 09:03 上传