pca学习笔记

主成分分析 \(\text{(PCA)}\) 算法学习笔记

\(\text{PCA}\) 应用场景

\(\text{deepseek}\) 给的回答

给定一个包含多个特征的高维数据集(例如学生的多科成绩、城市的各项经济指标等),希望利用 PCA 找到少数几个综合指标(主成分),既能最大程度保留原始数据的信息(方差),又能解释数据的主要结构,例如用于后续的降维、可视化或特征理解。

那么,开始!

算法与看法

将多个特征综合成少数几个新特征,仍能保证信息完整性的可行性与合理性,在于诸多特征是有相关性、受到现实因素制约的,如学生多个相近学科成绩一般相近。

于是可以科学地开始游戏。。。

首先获得我们要分析的数据集,记为矩阵 n 行 m 列 \(\mathbf{A}\),表示有 n 个样品,每个样品有 m 种特征因素,即 m 个维度的信息。

先做数据标准化消去量纲,并且去除不同特征数值波动范围差异的影响,毕竟我们不认为看到数值大的特征因素就认为这个维度很重要(因为可能数据在这个维度都比较大)。

\[{\mathbf{X}} = \frac{\mathbf{A} - \boldsymbol{\mu}}{\boldsymbol{\sigma}} \]

然后计算标准化后数据的协方差矩阵(明显找一手特征联系嘛)。

\[\mathbf{C}_{m \times m} = \frac{1}{m-1} \mathbf{X}_{n \times m}^{\top} \mathbf{X}_{n \times m} \]

证明很简单,用矩阵乘法展开验证即可,因为注意到标准化后均值为 0

这里的 \(m-1\) 有其统计学的意义,样本与总体联系,做的无偏估计,数学推导略(不是重点不要紧)

为什么要方差?因为我们去抓的主成分其实是找主导的、重要的因素

某维度数据方差大意味着这维度数据区分度大,也意味着这维度数据信息量大,是重要的

所以抓主成分就是找影响方差大的因素

这是理解层面

接下来从数学层面科学地去寻找

标准化后的数据集是 \(\mathbf{X}_{n \times m}\),假设我给出单位列向量 \(\mathbf{v}\),表示对这 \(m\) 维度的在意程度,那么 \(\mathbf{z} = \mathbf{X}_{n \times m}\mathbf{v}\) 就可以认为是每个样品的得分(总在意程度)

此时算方差 \(\operatorname{Var}(\mathbf{z})\),大的话意味着这些维度有很好的区分度,这些维度信息量大,是重要的

于是乎目标就是最大化方差,展开就可以证明下面这个式子

\[\operatorname{Var}(\mathbf{z})=\mathbf{v}^{\top} \mathbf{C} \mathbf{v} \]

\(\mathbf{v}\) 是特征向量时,带入特征分解 \(\mathbf{C} \mathbf{v} = \lambda \mathbf{v}\) 得到 \(\operatorname{Var}(\mathbf{z})=\lambda\) !!!

于是乎我们把特征分解和方差最大化联系起来了!

从几何上看,我们把 \(\mathbf{C}\) 特征分解重构得到特征基矩阵,做矩阵变换就可以视为特征基在其方向上的拉伸了,也就是特征分解的几何意义

联系到 \(\operatorname{Var}(\mathbf{z})\) 不难发现主导因素就被数学变换成了特征向量及其特征值,也就是旧坐标基底线性变换成的新坐标基底,新坐标基底可以按照对应的特征值排序确定重要性。

这里将线性变换的系数称为载荷矩阵 \(\mathbf{V}\) ,是由 \(\mathbf{v}\)组成的,各个特征值比特征值之和称为贡献率。那么就可以清楚地看到主导因素,也可以降维压缩了,因为受现实因素制约数据得到的特征值一般前几个贡献率较大,后面的可以视为噪声等不关键的因素,这样降维损失的信息就极少了

这里引入样品得分的概念,其实就是新坐标向量

\[\mathbf{Z} = \mathbf{X} \mathbf{V} \]

将每个样品丢进新坐标系并标出来的图称为散点图,在这个图上很容易看到团簇,从坐标分析出数据的属性特征

额外可以设计综合得分,以便比较样品间的优劣,不过这适用于新坐标所有维度正向或负向方向决定的好与坏统一的情况。

下面是鞭策 \(\text{deepseek}\) 许久得到的北太天元代码

实现了读入 data.csv 带有标题的样本数据,每行为一个样品,每列为一种特征

综合得分采用各个坐标维度按贡献率占比相加,并按正态分布转化为直观的百分制

输出文件为 pca_results.csv 包括贡献率数据、载荷矩阵、综合得分排名

同时也输出了一张图片包括贡献率直方图、降维后(按照保留 \(95\%\) 确定最终维数,但散点图视情况留2~3维以直观观察)的散点图,但需手动保存(没有调教出正确的自动保存代码。。。)

北太天元 pca_code

%% PCA主成分分析 - 正态分布百分制版
clear; clc;

%% 1. 读取数据
fprintf('===== PCA分析开始 =====\n');
fprintf('当前文件夹:%s\n', pwd);

% 读取CSV
fid = fopen('data.csv', 'r');

% 读特征名
第一行 = fgetl(fid);
特征名 = strsplit(第一行, ',');
特征名 = 特征名(2:end);

% 读数据
样本名 = {};
数据 = [];
行号 = 0;
while ~feof(fid)
    行 = fgetl(fid);
    if isempty(行), continue; end
    行号 = 行号 + 1;
    列 = strsplit(行, ',');
    样本名{行号,1} = 列{1};
    数据(行号,:) = str2double(列(2:end));
end
fclose(fid);

[n_samples, n_features] = size(数据);
fprintf('样本数:%d,特征数:%d\n', n_samples, n_features);

%% 2. PCA计算
X_std = (数据 - mean(数据)) ./ std(数据);
[V, D] = eig(cov(X_std));
[eigvals, idx] = sort(diag(D), 'descend');
V = V(:, idx);
scores = X_std * V;
贡献率 = eigvals / sum(eigvals) * 100;

%% 3. 确定保留的主成分个数
cum_贡献率 = cumsum(贡献率);
n_keep = find(cum_贡献率 >= 85, 1);
if isempty(n_keep)
    n_keep = n_features;
end
fprintf('保留前%d个主成分(累积贡献率%.1f%%)\n', n_keep, cum_贡献率(n_keep));

%% 4. 综合得分计算
权重 = 贡献率(1:n_keep) / 100;
综合得分_raw = scores(:, 1:n_keep) * 权重;

% 正态分布百分制转换
mu = mean(综合得分_raw);
sigma = std(综合得分_raw);
综合得分_z = (综合得分_raw - mu) / sigma;  % 标准化到N(0,1)

% 使用误差函数计算正态分布累积概率
综合得分_norm = 100 * 0.5 * (1 + erf(综合得分_z / sqrt(2)));

% 排序
[综合得分_sort, 排序索引] = sort(综合得分_norm, 'descend');
样本名_sort = 样本名(排序索引);

%% 5. 保存所有结果到一个CSV文件
fid_csv = fopen('pca_results.csv', 'w');

% 5.1 贡献率部分
fprintf(fid_csv, '主成分贡献率\n');
fprintf(fid_csv, '主成分,贡献率(%%),累积贡献率(%%)\n');
for i = 1:n_features
    fprintf(fid_csv, 'PC%d,%.2f,%.2f\n', i, 贡献率(i), cum_贡献率(i));
end
fprintf(fid_csv, '\n');

% 5.2 载荷矩阵部分
fprintf(fid_csv, '载荷矩阵\n');
fprintf(fid_csv, '特征名');
for i = 1:n_features
    fprintf(fid_csv, ',PC%d', i);
end
fprintf(fid_csv, '\n');
for i = 1:n_features
    fprintf(fid_csv, '%s', 特征名{i});
    for j = 1:n_features
        fprintf(fid_csv, ',%.4f', V(i,j));
    end
    fprintf(fid_csv, '\n');
end
fprintf(fid_csv, '\n');

% 5.3 样品综合得分排名部分(正态分布百分制)
fprintf(fid_csv, '样品得分排名\n');
fprintf(fid_csv, '排名,样本名,综合得分,百分制得分\n');
for i = 1:n_samples
    idx = 排序索引(i);
    fprintf(fid_csv, '%d,%s,%.4f,%.2f\n', ...
        i, 样本名{idx}, 综合得分_raw(idx), 综合得分_sort(i));
end

fclose(fid_csv);
fprintf('所有结果已保存到 pca_results.csv\n');

%% 7. 画图
figure('Position', [100, 100, 1200, 500]);

% 子图1:贡献率柱状图
subplot(1,2,1);
bar(贡献率, 'FaceColor', [0.3, 0.6, 0.9]);
xlabel('主成分编号');
ylabel('贡献率 (%)');
title('主成分贡献率');
grid on;
% 标出保留位置
hold on;
plot([n_keep+0.5, n_keep+0.5], [0, max(贡献率)], 'r--', 'LineWidth', 1.5);
text(n_keep+0.5, 2, sprintf('保留%d个', n_keep), ...
     'Color', 'red', 'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom');
% 标数值
for i = 1:n_features
    text(i, 贡献率(i)+1, sprintf('%.1f%%', 贡献率(i)), ...
         'HorizontalAlignment', 'center', 'FontSize', 9);
end

% 子图2:散点图(颜色用正态百分制)
subplot(1,2,2);
if n_keep >= 3
    scatter3(scores(:,1), scores(:,2), scores(:,3), 50, 综合得分_norm, 'filled');
    xlabel('PC1');
    ylabel('PC2');
    zlabel('PC3');
    title('样本得分分布(颜色:正态分布百分制)');
else
    scatter(scores(:,1), scores(:,2), 50, 综合得分_norm, 'filled');
    xlabel('PC1');
    ylabel('PC2');
    title('样本得分分布(颜色:正态分布百分制)');
end
colorbar;
grid on;

fprintf('\n===== PCA分析完成 =====\n');
fprintf('请手动保存图片(文件菜单 -> 另存为)\n');

实战演练

还没结束!!

经过 \(\text{deepseek}\) 老师的指点我找到了 AirQualityUCI.csv 作为 data.csv 来对代码进行测试

📊 UCI Air Quality Dataset 数据集档案

  • 数据名称:UCI Air Quality Dataset
  • 采样地点:意大利某城市,一个污染严重的道路旁
  • 采样时间:2004年3月到2005年2月(整整一年,每小时记录一次)
  • 样本数量:共 9358 条小时级别记录
  • 特征数量:除去日期时间,有 13 个可用于分析的数值特征
  • 缺失值标记:用 -200 代表数据缺失或无效

当然要修复一下以适配我上面的代码喽

依旧鞭笞 \(\text{deepseek}\)

北太天元 data_fix

%% 北太天元 - 修复空气质量数据格式
clear; clc;
fprintf('===== 开始修复数据格式 =====\n');

%% 1. 打开文件
fid = fopen('data.csv', 'r');
if fid == -1
    error('无法打开 data.csv,请确认文件在当前文件夹');
end

% 创建新文件
fid_new = fopen('data_fixed.csv', 'w');

%% 2. 读取并处理表头
% 读第一行
header_line = fgetl(fid);
% 按分号分割
header_parts = split(header_line, ';');
% 提取第3到15个作为特征名
feature_names = header_parts(3:15);

% 写入新表头
fprintf(fid_new, 'SampleName');
for i = 1:length(feature_names)
    fprintf(fid_new, ',%s', feature_names{i});
end
fprintf(fid_new, '\n');

%% 3. 处理数据行
line_count = 0;
empty_count = 0;
error_count = 0;

while ~feof(fid)
    % 读一行
    current_line = fgetl(fid);
    
    % 跳过空行
    if isempty(current_line) || strcmp(current_line(1), ';')
        empty_count = empty_count + 1;
        continue;
    end
    
    % 按分号分割
    parts = split(current_line, ';');
    
    % 检查是否有足够的数据列
    if length(parts) < 15
        error_count = error_count + 1;
        continue;
    end
    
    line_count = line_count + 1;
    
    % 合并日期和时间作为样本名
    date_str = parts{1};
    time_str = parts{2};
    % 把时间里的点换成冒号
    time_str = strrep(time_str, '.', ':');
    sample_name = [date_str, ' ', time_str];
    
    % 写入样本名
    fprintf(fid_new, '%s', sample_name);
    
    % 写入13个数值列
    for j = 3:15
        val_str = parts{j};
        if isempty(val_str)
            % 空值写NaN
            fprintf(fid_new, ',NaN');
        else
            % 把逗号换成点(欧洲数字格式)
            val_str = strrep(val_str, ',', '.');
            fprintf(fid_new, ',%s', val_str);
        end
    end
    fprintf(fid_new, '\n');
    
    % 显示进度
    if mod(line_count, 5000) == 0
        fprintf('已处理 %d 行...\n', line_count);
    end
end

%% 4. 关闭文件
fclose(fid);
fclose(fid_new);

%% 5. 显示统计信息
fprintf('\n===== 转换完成 =====\n');
fprintf('有效数据行数:%d\n', line_count);
fprintf('跳过的空行:%d\n', empty_count);
fprintf('格式错误行:%d\n', error_count);
fprintf('生成文件:data_fixed.csv\n');
fprintf('文件保存在:%s\n', pwd);

%% 6. 验证新文件
fprintf('\n===== 验证新文件 =====\n');
try
    % 尝试读取前几行看看
    data_test = readtable('data_fixed.csv', 'FileType', 'text', 'ReadVariableNames', true);
    fprintf('成功读取!\n');
    fprintf('样本数:%d\n', height(data_test));
    fprintf('变量数:%d\n', width(data_test));
    fprintf('变量名:\n');
    for i = 1:width(data_test)
        fprintf('  %d. %s\n', i, data_test.Properties.VariableNames{i});
    end
    fprintf('\n前3个样本:\n');
    disp(head(data_test, 3));
catch ME
    fprintf('验证失败:%s\n', ME.message);
end

好吧实测中我直接中断了 验证新文件 这个过程

AirQualityUCI_pca.png 如下

AirQualityUCI_pca

按照 \(95\%\) 的标准前四个维度足以,记为 \(\text{PC1,2,3,4}\)

表:主成分载荷矩阵

特征名 PC1 PC2 PC3 PC4
CO(GT) 0.0041 0.3619 0.5284 0.0609
PT08.S1(CO) -0.3638 0.1418 -0.1424 0.0541
NMHC(GT) -0.0383 0.0766 0.0442 0.9615
C6H6(GT) -0.3806 -0.0798 0.1047 -0.0531
PT08.S2(NMHC) -0.3378 0.1962 -0.2606 0.0073
NOx(GT) -0.0373 0.5279 0.1129 -0.1737
PT08.S3(NOx) -0.1431 -0.3805 0.5017 0.0641
NO2(GT) -0.0169 0.4828 0.3659 -0.0425
PT08.S4(NO2) -0.3374 0.0359 -0.2323 0.1214
PT08.S5(O3) -0.2986 0.2756 -0.2828 -0.0252
T -0.3610 -0.1549 0.1408 -0.0610
RH -0.3466 -0.1361 0.1766 -0.0798
AH -0.3628 -0.1479 0.1906 -0.0668

绝对值较大就意味着特征对新维度影响大

下面借助 \(\text{deepseek}\) 大力分析一手

主成分分析结果汇总

各主成分载荷特征

主成分 贡献率 主要载荷变量 正负含义 解释
PC1 50.65% PT08.S1、S2、S4、S5
T、RH、AH(均负载荷)
负值 → 污染严重
正值 → 空气清洁
夏季高温高湿型污染
PC2 22.65% 正:NOx、NO2、CO
负:PT08.S3
正值 → 新鲜交通污染
负值 → 二次转化
污染类型区分
PC3 10.84% CO、PT08.S3、NO2(均正) 正值 → 过程活跃
负值 → 过程不活跃
光化学反应程度
PC4 7.94% NMHC (0.96) 正值 → NMHC异常
负值 → NMHC正常
特殊排放事件

判断规则总结

主成分 污染方向 清洁方向 核心判断依据
PC1 ⬇️ 越负越污染 ⬆️ 正值 污染变量均为负载荷
PC2 ⬆️ 正值(一次)
⬇️ 负值(二次)
➡️ 接近0 正负代表污染类型
PC3 ⬆️ 正值 ⬇️ 负值 关键污染物载荷为正
PC4 ⬆️ 正值 ⬇️ 负值 NMHC 载荷 0.96

:PC1 的绝对值大小代表污染严重程度;PC2-PC4 的正负主要代表污染类型或过程状态。

于是在这个例子里综合得分没啥用

散点图聚集情况也很明显直观(虽然我不太会读

总之无论如何算成功了

完结撒花!!!

posted @ 2026-03-16 23:22  leiyuanze  阅读(0)  评论(0)    收藏  举报