主成分分析（数据分析课本例3.9.1）

Python代码：

# -*- coding: utf-8 -*-
"""
Created on Tue Feb  6 14:37:21 2018
E-mail = Eric2014_Lv@sjtu.edu.cn
@author: DidiLv
"""
# 项目内容：线性统计模型（线性回归与方差分析）例3.9.1
import numpy as np
from sklearn import preprocessing
from numpy import linalg as la
# 导入数据


# 导入x
X1 = np.array([ [149.3, 4.2, 108.1],
                [ 161.2, 4.1, 114.8],
                [171.5, 3.1, 123.2],
                [175.5, 3.1, 126.9],
                [180.8, 1.1, 132.1],
                [190.7, 2.2, 137.7],
                [202.1, 2.1, 146.0],
                [212.4, 5.6, 154.1],
                [226.1, 5.0, 162.3],
                [231.9, 5.1, 164.3],
                [239.0, 0.7, 167.6]])
# 导入y
Y1 = np.array([15.9, 16.4, 19.0, 19.1, 18.8, 20.4, 22.7, 26.5, 28.1, 27.6, 26.3])


# 中心化标准化直接用相关的包进行处理
X = preprocessing.scale(X1)
Y = preprocessing.scale(Y1)
# 特征值分解
U, Sigma, VT = la.svd(np.dot(X.T, X))
#将第三个特征值为0
Sigma[2] = 0.0
#Z就是主成分
Z = np.dot(X,U)
#我们选取前两个主成分：就是前两列
Z1 = Z[:,:2]
Z1TZ1_inv = la.inv(np.dot(Z1.T, Z1))
# 这是对 alpha1的估计
alpha1 = np.dot(np.dot(Z1TZ1_inv,Z1.T),Y)
# 下面就是给alpha1某尾加上0
alpha = np.append(alpha1,[0.0])
#这是最终的beta的数值
beta = np.dot(U,alpha)
'''
# 因为到现在都是中心化和标准化的数据进行的回归，
# 还有一步就是利用69页下面的那个经验回归方程这样的形式还原成3.9.1的例子的形式，
# 你完成吧！
'''

Matlab代码：

clear
%%  导入数据（未经过中心化标准化的数据）
% 导入x
X1 = [ 149.3, 4.2, 108.1;
      161.2, 4.1, 114.8;
      171.5, 3.1, 123.2;
      175.5, 3.1, 126.9;
      180.8, 1.1, 132.1;
      190.7, 2.2, 137.7;
      202.1, 2.1, 146.0;
      212.4, 5.6, 154.1;
      226.1, 5.0, 162.3;
      231.9, 5.1, 164.3;
      239.0, 0.7, 167.6;];
  % 导入y
 Y1 = [15.9, 16.4, 19.0, 19.1, 18.8, 20.4, 22.7, 26.5, 28.1, 27.6, 26.3]';
 %% 中心化标准化： zscore函数
[m,n]=size(X1);  
X = zscore(X1); 
Y = zscore(Y1); 
%% 特征值分解： Matlab里面的svd对于对称矩阵的时候就是特征值分解
[U, Sigma, V] = svd(X'*X); % 分解后的结果： X'X = U*Sigma*V' ; 
Sigma(3,3)=0; %将第三个特征值为0
Z = X*U; % Z就是主成分
Z1 = Z(:,1:2); % 我们选取前两个主成分：就是前两列
Z2 = Z(:,3); % 这是第三个主成分（第三列），由于对应的特征值为0，所以舍掉，在后面没有用到，只是写到这里符合课本的顺序。
alpha0 = mean(Y); % 这是对alpha0的估计
alpha1 = inv(Z1'*Z1)*Z1'*Y % 这是对 alpha1的估计
alpha = [alpha1; 0] % 我们看到alpha的第三个量为零，实际上就是把第三个主成分对应的系数改为了0
beta = U*alpha %  
%% 还有一步就是利用69页下面的那个经验回归方程这样的形式还原成3.9.1的例子的形式，你完成吧！
%