网站建设公司推荐金石下拉网络,双语外贸网站源码,游戏制作流程,wordpress发布视频文章目录Part.I 预备知识Chap.I 一些概念Chap.II 主成分分析Chap.III Matlab 函数 randnChap.IV Matlab 函数 pcaPart.II 置信椭圆的含义Chap.I 一个 Matlab 实例Sec.I 两个不相关变量的特征Sec.II 两个相关变量的特征Chap.II 变换阵 (解相关矩阵) 的求解ReferencePart.I 预备知…
文章目录Part.I 预备知识Chap.I 一些概念Chap.II 主成分分析Chap.III Matlab 函数 randnChap.IV Matlab 函数 pcaPart.II 置信椭圆的含义Chap.I 一个 Matlab 实例Sec.I 两个不相关变量的特征Sec.II 两个相关变量的特征Chap.II 变换阵 (解相关矩阵) 的求解ReferencePart.I 预备知识
Chap.I 一些概念
首先要了解一下下面的概念
点估计设总体X的分布函数的形式已知但它的一个或多个参数未知借助于总体X的一个样本来估计总体未知参数的值的问题称为参数的点估计问题。对于一个未知量人们在测量或计算时常不以得到近似值为满足还需估计误差即要求知道近似值的精确程度(亦即所求真值所在的范围)。类似地对于未知参数 θ\thetaθ除了求出它的点估计 θ^\hat \thetaθ^ 外我们还希望估计出一个范围并希望知道这个范围包含参数 θ\thetaθ 真值的可信程度。这样的范围通常以区间的形式给出同时还给出此区间包含参数 θ\thetaθ 真值的可信程度。这种形式的估计称为区间估计这样的区间即所谓置信区间。置信区间Confidence interval 设总体 X 的分布函数 F(x;θ)F(x;\theta)F(x;θ) 含有一个未知参数 θ\thetaθθ∈Θ\theta\in\Thetaθ∈Θ Θ\ThetaΘ 是 θ\thetaθ 可能取值的范围对于给定值 α\alphaα0α10\alpha10α1若由来自 X 的样本 X1,X2,⋯,XnX_1,X_2,\cdots,X_nX1,X2,⋯,Xn 确定的两个统计量 θ‾θ‾(X1,X2,⋯,Xn)\underline \theta\underline \theta(X_1,X_2,\cdots,X_n)θθ(X1,X2,⋯,Xn) 和 θ‾θ‾(X1,X2,⋯,Xn)(θ‾θ‾)\overline \theta\overline \theta(X_1,X_2,\cdots,X_n)\ (\underline \theta\overline \theta)θθ(X1,X2,⋯,Xn) (θθ) 对于任意 θΘ\theta\ThetaθΘ 满足 P{θ‾(X1,X2,⋯,Xn)θθ‾(X1,X2,⋯,Xn)}≥1−αP\{\underline \theta(X_1,X_2,\cdots,X_n)\theta\overline \theta(X_1,X_2,\cdots,X_n)\}\ge1-\alphaP{θ(X1,X2,⋯,Xn)θθ(X1,X2,⋯,Xn)}≥1−α则称随机区间 (θ‾,θ‾)(\underline \theta,\overline \theta)(θ,θ) 是 θ\thetaθ 的置信水平为 1−α1-\alpha1−α 的置信区间θ‾\underline \thetaθ 和 θ‾\overline \thetaθ 分别称为置信水平为 1−α1-\alpha1−α 的双侧置信区间的置信下限和置信上限1−α1-\alpha1−α 称为置信水平。上面给出了置信区间严格的数学定义通俗的讲置信区间展现的是这个参数的真实值有一定的带来落在测量结果的程度其给出的是被量测参数的量测值的可信程度置信水平一般常用 95%95\%95% 置信水平。68−95−99.768-95-99.768−95−99.7 法则或者叫做 3σ3\sigma3σ 原则。对于正态分布68%68\%68% 的数据分布在距离均值1个方差的范围95%95\%95% 的数据分布在距离均值2个方差的范围99.7%99.7\%99.7% 的数据分布在距离均值3个方差的范围内。
Chap.II 主成分分析
主成分分析 (Principal Component Analysis, PCA) 是一种统计方法通过正交变换将一组可能存在相关性的变量转换为一组线性不相关的变量转换后的这组变量叫主成分。PCA 的主要应用有降维、特征提取、去噪、故障检测。
PCA 的主要思想是将 nnn 维特征映射到 kkk 维上这 kkk 维是全新的正交特征也被称为主成分是在原有 nnn 维特征的基础上重新构造出来的 kkk 维特征。PCA 的工作就是从原始的空间中顺序地找一组相互正交的坐标轴新的坐标轴的选择与数据本身是密切相关的。其中第一个新坐标轴选择是原始数据中方差最大的方向第二个新坐标轴选取是与第一个坐标轴正交的平面中使得方差最大的第三个轴是与第1,2个轴正交的平面中方差最大的。依次类推可以得到 nnn 个这样的坐标轴。通过这种方式获得的新的坐标轴我们发现大部分方差都包含在前面 kkk 个坐标轴中后面的坐标轴所含的方差几乎为0。于是我们可以忽略余下的坐标轴只保留前面 kkk 个含有绝大部分方差的坐标轴。事实上这相当于只保留包含绝大部分方差的维度特征而忽略包含方差几乎为0的特征维度实现对数据特征的降维处理。
Chap.III Matlab 函数 randn
下面的实例中会用到 Matlab 中生成正态分布随机数的函数randn()下面对其用法做简要的介绍。
randn产生均值为0、方差为1的高斯白噪声常用的用法有下面几种randn(n)产生一个 n*n 的随机数方阵randn(m,n)产生一个 m*n 的随机数矩阵 除此之外还有Matlab 关于随机数的生成函数还有
rand产生均值为0.5、幅度在0~1之间的伪随机数。randperm(n)产生1到n的均匀分布随机序列。normrnd(a,b,c,d)产生均值为a、方差为b大小为c*d服从正态分布的随机矩阵。 下面对randn做一个简单的验证产生两列数据分别画出它们的概率分布看看图形是否符合正态分布。代码和绘制结果如下可以看出结果符合预期。
%% 检验 randn 函数
function rTest()
datarandn(600,2);
figure
hold on
subplot(121)
t2;
pcapaplot(data(:,1),[-t,t]);
subplot(122)
pcapaplot(data(:,2),[-t,t]);
endChap.IV Matlab 函数 pca
function [coeff, score, latent, tsquared, explained, mu] pca(x,varargin)输入x 是样本N*PN个样本量P维 输出
COEFF返回N×P数据矩阵X的主成分系数。X的行对应于观测值列对应于变量。每列系数包含一个主成分的系数。各列按主成分方差latent降序排列。默认情况下PCA将数据居中并使用奇异值分解算法。对于非默认选项请使用名称/值对参数。SCORE返回主成分得分它是X在主成分空间中的表示。score的行对应观察值列对应主成分。中心数据可以用SCORE*COEFF重建。LATENT返回每个主成分方差即X的协方差矩阵的特征值特征值从大到小进行排序。TSQUARED返回X中每个观测值的Hotelling T平方统计值。PCA使用所有主分量计算TSQUARED在整个空间中计算即使请求的组件较少请参见下面的“NumComponents”选项。对于缩小空间中的TSQUARED使用MAHALSCORESCORE。EXPLAINED返回一个向量其中包含每个主成分方差占总方差的百分比。MUCentered设置为true时返回估计的平均值MU设置为false时返回所有零。
Part.II 置信椭圆的含义
对于一维的情况我们可以用置信区间来描述一个值落在某个区间区域的概率置信水平这个区间用几何表示是一个线段对于二维的情况几何意义上就从一条直线变成了一个平面类似于置信区间便有了置信椭圆confidence ellipse的概念。 下面是笔者的思考一些自说自话罢了。 置信椭圆的概念是从置信区间类比而来的严格意义上讲它描述的仅仅是二维的情况对于三维可以称之为置信椭球对于更高维可以称之为置信超椭球自己瞎起的~但是很多情况下对于二维或者更高维很多学者也将其称为置信椭圆。为什么叫『椭圆』不叫『圆』『长方形』呢因为点的分布从图形上看可以用一个椭圆将其包裹起来所以叫椭圆。『误差椭圆』这个概念又是怎么来的对于二维情况我们进行多次量测的目的是为了确定一个二维参数的取值但是我们又不能 100%100\%100% 确定我们给出的值就是这个二位参数的真值但是我们一般可以说这个二维参数有多少概率会落在某个范围内。在给定概率的情况下这个范围越小说明我们对我们的量测越有信心我们量测的精度越高误差越小。这个范围可以说是我们给出参数的精度或误差所以便有了『误差椭圆』这个概念这么说『精度椭圆』这个名词是不是也顺利成章了 Chap.I 一个 Matlab 实例
其实写这篇博文的初衷是为了搞清楚『为什么两个随机变量不相关其误差椭圆的长轴平行于坐标轴』。还是看个例子用数据说话吧。在看下面的内容之前首先要明确下面几点
我们通过 randn 函数来生成了两列均值为 0 方差为 1 的随机数我们对这列随机数的理解为每一列数据都表示对一个『代求量』的多次量测显然我们已经假设我们的量测存在一个服从高斯分布的噪声并且我们的『代求量』真值是0。我们也可以将每一列数据理解为对一个『随机变量』的多次抽样或观测。randn 函数虽然一次生成了两列数据但是这两列数据是毫无关联的当然也可以一次生成一列数据分两次生成换言之这两列数据代表的两个『随机变量』是不相关的。
Sec.I 两个不相关变量的特征
首先用randn生成了两列均值为0、方差为1的随机数据数据量是300*n为了探究不同数据量的区别因此画了6幅图如下数据量分别是300*nn1,2,3,4,5,6。图中len表示数据量的大小a,b 分别表示椭圆长轴和短轴长度红色和蓝色箭头好吧红色可能看不到箭头分别表示长轴和短轴方向。可以看到画出来的图形几近与圆数据量越多它越接近于圆。但是它的轴并不平行于坐标轴可以这样解释因为这两列数据的方差都是1所以长短轴理论上应该相等所以画出来的图形理论上应该是个圆实际情况也是数据量越多它越接近于圆所以其指向就无所谓了。 为了让两列数据的方差不同我们将第一列数据乘上了2不同数据量的绘图结果如下。可以看到它确实是长短轴平行于坐标轴的椭圆 涉及的绘图代码如下
%% 两个不相关变量的特性
function Test1(wid,hei,a,b)
figure
pmag50; %the size of margin pix
marx(1:b1,1)3;marx(1)8;
mary(1:a1,1)3;mary(1)8;
sumxsum(marx);sumysum(mary);
marxmarx./sumx;marymary./sumy;
lefypmag/hei/mary(1);
lefxpmag/wid/marx(1);
marxmarx.*lefx;
marymary.*lefy;
wid1(1-lefx)/b;
hei1(1-lefy)/a;xt-4:2:4;xtl1sprintfc(%g,xt);xtl2kcell(size(xtl1));xl[-5 5];
yt1-4:2:4;ytl1sprintfc(%g,yt1);ytl2kcell(size(ytl1));yl[-5 5];set(gcf,position,[0 0 wid hei])
for i1:amagy(a-i)*hei1sum(mary(1:end-i));for j1:bmagx(j-1)*wid1sum(marx(1:j));n(i-1)*bj;subplot(a,b,n);set(gca,position,[magx magy wid1 hei1],box,on);% plot here!len300;datarandn(n*len, 2); % get the data PCA(data);set(gca,XLim,xl,XTick,xt,XTicklabel,xtl2);set(gca,YLim,yl,YTick,yt1,YTicklabel,ytl2);if j1ylabel(Ylab);set(gca,YTick,yt1,YTicklabel,ytl1);endif iaxlabel(Xlab);set(gca,XTick,xt,XTicklabel,xtl1);endend
end
end% get a null cell which is a*b dims
function datakcell(m)
am(1);bm(2);
datacell(a,b);
for i1:afor j1:bdata(i,j)cellstr(num2str(data{i,j}));end
end
end%% PCA
function PCA(data)
lensize(data,1);
data(:,1)2*data(:,1);
center mean(data);
[coeff, ~, latent, ~, ~] pca(data);% r1 r2 为自定义的向量大小参数(模)
r1 6;
r2 3;
% p1 p2 为第一主轴和第二主轴上的点
p1 r1*coeff(:, 1)center;
p2 r2*coeff(:, 2)center;% 主轴方向与X轴之间的夹角
angle cart2pol(coeff(1, :), coeff(2, :))*180/pi;
beta angle(1, 1);
% 置信椭圆坐标以 95% 为例
semimajor sqrt(latent(1, 1)); % 长轴长度一半
semiminor sqrt(latent(2, 1)); % 短轴长度一半
alpha linspace(0, 360, 2000);
% 卡方分布表
% https://people.richland.edu/james/lecture/m170/tbl-chi.html
% level 4.605; % 90%
level 5.991; % 95%
% level 9.210; % 99%
% 椭圆坐标点
ellipse_X center(1, 1)sqrt(level)*(semimajor*cosd(alpha)*cosd(beta)-...semiminor*sind(alpha)*sind(beta));
ellipse_Y center(1, 2)sqrt(level)*(semimajor*cosd(alpha)*sind(beta)...semiminor*sind(alpha)*cosd(beta));%% 可视化
% figure
hold on
box on
grid on
% 原始数据
scatter(data(:, 1), data(:, 2), 15, LineWidth, 1.2,...MarkerEdgeColor, [151, 138, 189]/255,...MarkerFaceColor, [151, 138, 189]/255);
xlim([-5, 5]);
ylim([-5, 5]);
set(gca, linewidth, 1.5)% 置信椭圆
plot(ellipse_X, ellipse_Y, Color, [0, 102, 255]/255,...LineStyle, -, LineWidth, 3),% 第一主轴方向
arrow_1 annotation(arrow, Color, [162, 20, 47]/255,...HeadStyle, cback2, LineWidth, 3, HeadWidth, 20, HeadLength, 20);
arrow_1.Parent gca;
arrow_1.X [center(1, 1), p1(1, 1)];
arrow_1.Y [center(1, 2), p1(1, 2)]; % 第二主轴方向
arrow_2 annotation(arrow, Color, [0, 114, 189]/255,...HeadStyle, cback2, LineWidth, 3, HeadWidth, 20, HeadLength, 20);
arrow_2.Parent gca;
arrow_2.X [center(1, 1), p2(1, 1)];
arrow_2.Y [center(1, 2), p2(1, 2)]; % 中心点
plot(center(1, 1), center(1, 2),...Marker, o,...MarkerSize, 8,...MarkerEdgeColor, [162, 20, 47]/255,...MarkerFaceColor, [162, 20, 47]/255);% title(主轴方向和置信椭圆, FontSize, 16, FontWeight, bold)
strsprintf(len%d\na%.2f\nb%.2f,len,semimajor,semiminor);
text(2, -3, str, FontSize, 16, FontWeight, bold)
axis equal
endSec.II 两个相关变量的特征
为了让两个随机变量相关我们记randn生成的两列数据为 X[a,b]X[a,b]X[a,b]其中 aaa 的均值为0方差为2bbb 的均值为0方差为1将其乘上一个矩阵得到新的数据 YYY Y[a′,b′][a,b]⋅[11−12]X⋅ZY[a,b][a,b]\cdot \left[ \begin{array}{ccc} 1 1 \\ -1 2 \\ \end{array} \right]X\cdot ZY[a′,b′][a,b]⋅[1−112]X⋅Z 也就是说a′a−b,b′a2baa-b,\ ba2ba′a−b, b′a2b然后我们绘制其置信椭圆如下图所示 其代码如下
%% 两个相关变量的特征
function Test2()
len300;
datarandn(5*len, 2); % get the data
data(:,1)2*data(:,1); % 第一列方差变为2
data1zeros(size(data));
Z[1 1;-1 2];
data1data*Z;
% data1(:,1)data(:,1)-data(:,2);
% data1(:,2)data(:,1)2*data(:,2);
PCA(data1);
end我们能让两个不相关的随机变量通过乘以一个矩阵变成相关的同样的我们也可以通过乘一个矩阵使它们『解相关』。对于这组数据我们已经知道了使它俩变成相关的矩阵 ZZZ那么再乘 Z−1Z^{-1}Z−1 显然可以解相关但是对于一个我们不知道变换矩阵对应此例也就是上文的 Z−1Z^{-1}Z−1 的数据我们一般用下面的方法来求『变换阵』。 PS之前的我还以为用 LU 分解来求变换阵并且认为求解变换阵的方法是『高斯变换』变换阵被称为『高斯变换阵』。但是并不是是通过特征值和特征向量来求的。 Chap.II 变换阵 (解相关矩阵) 的求解
现在的问题是已知两个随机变量的许多量测两列数据而且这两个随机变量是相关的怎么找到一个矩阵使得这两列数据乘上这个矩阵之后这两列数据就不相关了。笔者觉得有两种方法一种是几何意义上的将置信椭圆旋转一下使得其长短轴平行于坐标轴旋转在代数意义上乘以一个矩阵另一种是纯代数通过矩阵变换找到这样的旋转矩阵。
实际上上述问题用主成分分析找主元、降相关就可以解决。 基于特征值分解协方差矩阵实现 PCA 算法步骤
输入数据集 X{x1,x2,⋯,xn}X\{x_1,x_2,\cdots,x_n\}X{x1,x2,⋯,xn}需要将其降到 kkk 维去中心化每一维特征减去各自的均值计算方差-协方差矩阵 A1nXTXA\frac{1}{n}X^TXAn1XTX 这里分母是 nnn 或 n−1n-1n−1 对求出的特征向量没有影响求出 AAA 的特征值和特征向量对特征值从大到小排序选择其中最大的 kkk 个。然后将其对应的 kkk 个特征向量分别作为行向量组成特征向量矩阵 PPP将数据转换到 kkk 个特征向量构建的新空间中即 YPXYPXYPX 下面对此种方法求『解相关矩阵』进行一个简单的证明。
首先要明确我们的目的是将 XXX 的协方差矩阵 AAA 转换为一个对角阵即 XTXAVTΛVX^TXAV^T\Lambda VXTXAVTΛV也就是令 (XV−1)T(XV−1)Λ(XV^{-1})^T(XV^{-1})\Lambda(XV−1)T(XV−1)Λ即寻求一个转换阵 ZV−1ZV^{-1}ZV−1使得 BTBΛ(BXZ)B^TB\Lambda\ (BXZ)BTBΛ (BXZ)关于矩阵分解的方法中哪种分解方法可以使得矩阵的分解结果中含有一个对角阵呢奇异值SVD和特征分解分解可以用这两种分解中的哪种都可以特征分解AVΛV−1AV\Lambda V^{-1}AVΛV−1 并且协方差阵是一个实对称阵对于实对称阵其不同特征值对应的特征向量正交。又因为特征向量伸缩乘一个常数之后还是特征向量所以我们很容易构造出来一个特征向量矩阵 VVV使得 VTV−1V^TV^{-1}VTV−1.基于上述有 XTXVΛV−1X^TXV\Lambda V^{-1}XTXVΛV−1. 所以我们取转换阵 ZVZVZV 令 YXZYXZYXZYTYVTXTXVVTVΛV−1VV−1VΛV−1VΛY^TYV^TX^TXVV^TV\Lambda V^{-1}VV^{-1}V\Lambda V^{-1}V\LambdaYTYVTXTXVVTVΛV−1VV−1VΛV−1VΛ。这样利用转换阵 VVV 便将协方差矩阵变成了对角阵即实现了『解相关』。 基于上述步骤对于相关变量得到了其变换阵 X1X1X1绘图结果如下左上是原始数据左下和右上对应中可行的转换矩阵右下啥也不是。我们发现得到的变换阵 X1X1X1 并不是 Z−1Z^{-1}Z−1 DataYDataX⋅ZData_YData_X\cdot ZDataYDataX⋅Z这说明变换阵不唯一 代码如下
%% PCA 求 变换矩阵 并绘图
function Test3()
len300;
datarandn(5*len, 2); % get the data
data(:,1)2*data(:,1); % 第一列方差变为2
Z[1 1;-1 3];
data1data*Z;
Adata1*data1/(5*len);
[X,D]eig(A); % 求特征值 D 和特征向量 X; X 各列是相应的特征向量
% A*X(:,1)-D(1,1)*X(:,1)
data2data1*X;
X1[X(:,2),X(:,1)]; % 因为第一个特征值比较大所以优先将第一个特征向量放在前面
data3data1*X1;
data4data1*X1;
% Plot
subplot(2,2,1);
PCA(data1); % 原始数据
subplot(2,2,2);
PCA(data2); % 去相关之后的 y 轴长
subplot(2,2,3);
PCA(data3); % 去相关之后的 x 轴长
subplot(2,2,4);
PCA(data4); % 啥也不是
endPS实际上笔者使用的 PCA 函数中已经做了主成分分析求出了样本的方差-协方差矩阵里面用的是Matlab自带的pca函数为了更详细地了解 PCA 的过程笔者又做了点额外工作PCA 函数在本博文中仅仅发挥一个画图的功能。 Reference
主成分分析 (PCA) 的主轴和置信椭圆可视化主成分分析PCA原理详解PCA的数学原理置信椭圆与R画法Matlab 让多图排版更美观矩阵的各种分解汇总注上面的代码已经很全了笔者也提供了测试时所使用的.m文件戳我下载。
