MATLAB数据分析方法第2章数据描述性分析

2.1基本统计量与数据可视化

1.均值、中位数、分位数、三均值

均值、中位数：mean(A)、media(A)

分位数：prctile(A,P)，P∈[0,100]

prctile(A,[25,50,75]) %求A的下、中、上分位数

三均值：

w=[0.25,0.5,0.75];

SM=w*prctile(A,[25,50,75])

%例：计算安徽16省市森林资源统计量
A=xlsread('senlin.xls','sheet1')
M=mean(A); %均值, 
MD=median(A); %中位数
SM=[0.25,0.5,0.25]*prctile(A,[25,50,75]); %三均值
[M;MD;SM]

2.方差、标准误、变异系数

方差：var(A,flag),flag默认0表示修正的方差，取1为未修正

标准差：std(A,flag),同上

变异系数：v=std(A)./abs(mean(A))

k阶原点矩、中心距：

ak=mean(A.^k)

bk=mean((A-mean(A)).^k)

%中心距系统命令bk=moment(A,k)

3.极差、四分位极差(上、下分位数之差)

R=rangr(A)

R1=iqr(A)

4.异常点判别（截断点）

XJ=parctile(A,[25])-1.5*R1

SJ=parctile(A,[75])+1.5*R1

5.偏度、峰度

偏度：sk=skewness(A,flag),默认1，取0为样本数据修正的偏度

峰度：ku=kurtosis(A,flg)-3,同上

2.1.2多维样本数据

协方差：cov(A)

相关系数：corr(A)

标准化：zscore(A)

2.1.3样本数据可视化

1.条形图

bar(x)%样本数据x的条形图，横坐标为1：length(x)

bar(x,y)%先把x和y一一对应，然后将x从小到大排序画图

2.直方图

hist(x,n)%数据x的直方图，n为组数，确省时n=10

[h,stats]=cdfplot(x)%x的经验分布函数图，stats给出数据最大最小值、中位数、均值、标准差

直方图基础上附加正态密度曲线

histfit(x)

histfit(x,nbins)%nbins指定bar个数，缺省时为x中数据个数的平方根

3.盒图，五个数值点组成：最小值、下四分位数、中位数、上四分位数、最大值。中间的盒子从Q1延申到Q3，盒子的直线标出中位数位置，盒子两端有直线往外延伸到最小数与最大数

boxplot(x)%矩阵x的每一列的盒图和“须”图，

4.阶梯图

stairs(x)%x的阶梯图

5.火柴棒图

stem(x)%离散数据的火柴棒图

%例：随机150个服从标准正态分布随机数，做出柱形图、直方图、阶梯图、火柴棒图
x = random('normal',0,1,[1,150]);   %产生服从标准正态分布随机数150个
subplot(2,3,1),bar(x),title('柱形图')
subplot(2,3,2),hist(x),title('直方图')
subplot(2,3,3),stairs(x),title('阶梯图')
subplot(2,3,4),stem(x),title('火柴棒图')
subplot(2,3,5),histfit(x),title('附正态密度曲线')
subplot(2,3,6),boxplot(x),title('盒图')

例：（X,Y）服从二维正态分布N(2,1;3,3;sqrt(3)/2),生成100对数据做出三点图画出二维正态分布的密度函数，绘制密度曲面图形
mu = [2 3];                  %输入均值向量
sa = [1 1.5; 1.5 3];          %输入协方差矩阵
Nr = mvnrnd(mu,sa,100);       %随机生成n=100的样本数据
scatter(Nr(:,1),Nr(:,2),'*'); %作样本数据平面散点图
%绘制密度曲面
figure(2)
v=sqrt(3)/2;                 %输入相关系数
x=-1:0.05:5;                 %横坐标的取值向量
y=-2:0.05:8;                 %纵坐标的取值向量
[X,Y]=meshgrid(x,y);         %生成网格点
T=((X-mu(1)).^2/sa(1,1)-2*v/sqrt(sa(1,1)*sa(2,2))*(X-mu(1)).*(Y-mu(2))+(Y-mu(2)).^2/sa(2,2));
Z=1/(2*pi)/sqrt(det(sa))*exp(-1/2/(1-3/4)*T);   %计算密度函数值
mesh(X,Y,Z)                  %绘制密度曲面图形

 

3.正态概率图和Q-Q图

(1)正态概率图：normplot(x)%样本数据对应图中的+，若+大都集中在红色参考线上则说明服从正态分布

若总体是非正态分布，可以类似绘制概率图（威布尔分布概率图）

weibplot(x)%若样本数据点基本散布在一条直线上，则数据服从该分布

(2)Q-Q图：

PD=fitdist(x,distname)%distname为Beta、Binomial、Exponential、Normal、Weibull

qqplot(x,PD)% PD省略时为标准正态分布

%例：自由度为8的卡方分布数据300个，绘制其正态概率图与卡方分布的Q-Q图
clear
s=rng;rng(s);
c1=chi2rnd(8,[300,1]); c2=sort(c1);      %模拟生成卡方分布样本
plot(c2,chi2pdf(c2,8), '+-');             %绘制卡方分布的密度曲线
title('卡方分布的密度曲线') ;legend('自由度n=8');
grid on									%打开网格
figure
pd=makedist('Gamma','a',4,'b',0.5)   %创建参数a=4,b=0.5的伽马分布
%pd=gamrnd(4,0.5,[300,1]);        
subplot(1,2,1),normplot(c1);        %绘制样本的正态概率图
subplot(1,2,2),qqplot(c1,pd);        %按指定分布绘制样本的Q-Q图
grid on

左图可看出是左偏，中图偏离参考线故不服从正态分布，右图贴合直线故服从伽马分布（自由度为8的卡方分布就是参数a=n/2,b=0.5的伽马分布）

2.2数据分布及其检验

1.经验分布函数

[h,stats]=cdfplot(x)%x的经验分布函数图，h表示曲线的环柄，stats给出数据最大值、最小值、中位数、均值、标准差

%例：生成服从标准正态分布的50个样本点，画出样本的经验分布函数图，并与理论分布函数比较
clear
X=normrnd (0,1,50,1);           %生成服从标准正态分布的50行1列样本点
[h,stats]=cdfplot(X);           %样本的经验分布函数图
hold on							%保持上图，后续图层叠加在上图
plot(-3:0.01:3, normcdf(-3:0.01:3,0,1), 'r')   %理论分布函数图
legend('样本经验分布函数Fn(x)', '理论分布函数Φ(x)','Location','NorthWest')

2.总体分布正态性检验

(1)JB检验：显著性水平alpha默认0.05

%以下两种方法均可，输出的H=0无法拒绝X服从正态分布；H=1拒绝。P为接受假设的概率值，P<alpha，则拒绝正态分布原假设。JBSTAT为测试统计量，CV是拒绝临界值，JBSTAT>CV拒绝。

H=jbtest(x,alpha)

[H,PJBSTAT,CV]=jbtest(x,alpha)

(2)KS检验：

h=kstest(x)

h= kstest(x,cdf)

[h,p,ksstat,cv]=kstest(x,cdf,alpha)

%h=0无法拒绝X服从正态分布

(3)Lilliefors检验(改进的KS检验)：

H=lillietest(x,alpha)

[H,P,LSTAT,CV]=lillietest(x,alpha)

%判断方法同（1）

%例：检验“中国银行”的股票的收盘价是否服从正态分布。
clear
a=xlsread('yhgspj.xls');  %读取收盘价数据
h1=jbtest(a(:,1))    %JB检验
h2=kstest(a(:,1))    %KS检验
h3=lillietest(a(:,1))   %改进KS检验
qqplot(a(:,1))

h1=1,h2=1,h3=1,三种检验均拒绝正态分布原假设，Q-Q图显示偏离参考线，不服从正态分布

2.2.2多维数据的正态分布检验

%例：三组数据有4项指标，三组视为一个总体，是否服从四维正态分布
clear
A=xlsread('jibing.xls','sheet1')
X=[A(:,1:4);A(:,5:8);A(:,9:12)];
[N,p]=size(X);
d=mahal(X,X);                % 计算马氏距离
d1=sort(d);                    % 从小到大排序
pt=[[1:N]-0.5]/N;               % 计算分位数
x2=chi2inv(pt,p);               % 计算X2t
plot(d1,x2','*',[0:12],[0:12],'-r')  % 作图

数据点基本落在直线上，不能拒绝正态分布原假设

3.多个总体协方差矩阵相等性检验

(1)两总体

%例：蠓虫分类
clear
apf=[1.14,1.78; 1.18,1.96;1.20,1.86;1.26,2.;1.28,2;1.30,1.96];
af=[1.24,1.72;1.36,1.74;1.38,1.64;1.38,1.82;1.38,1.90;1.40,1.70;1.48,1.82;1.54,1.82;1.56,2.08];
n1=6;n2=9;p=2;
s1=cov(apf);s2=cov(af);
s=((n1-1)*s1+(n2-1)*s2)/(n1+n2-2);             % 计算混合样本方差
Q1=(n1-1)*(log(det(s))-log(det(s1))-p+trace(inv(s)*s1));
Q2=(n2-1)*(log(det(s))-log(det(s2))-p+trace(inv(s)*s2));% 计算检验统计量观测值
ch2inv(0.95,3)                                %计算X2(3)临界值

Q1=2.5784<7.8145,Q2=0.7418<7.8145,无法拒绝两类总体协差阵相同的原假设

(2)多总体

%例：上上个例子三组数据的协方差矩阵是否相等？...
clear,clc
A=xlsread('jibing.xls','sheet1');      %输入样本数据
G1=A(:,1:4);                           %提取总体1的样本
G2=A(:,5:8);        
G3=A(:,9:12);
n=size(G1,1)+ size(G2,1)+ size(G2,1);  %计算总的样本容量
[n1,p]= size(G1);
k=3;
f=p*(p+1)*(k-1)/2;                      %统计量自由度
d=(2*p^2+3*p-1)*(k+1)/(6*(p+1)*(n-k));
s1=cov(G1);                             %协方差矩阵
s2=cov(G2);                             %协方差矩阵
s3=cov(G3);                             %协方差矩阵
s=(n1-1)*(s1+s2+s3)/(n-k);              %总体协方差估计
M=(n-k)*log(det(s))-19*(log(det(s1))+log(det(s2))+log(det(s3))); 
T=(1-d)*M
P0=1-chi2cdf(T,f)                      % 卡方分布概率

P=0.4374>0.1故无法拒绝三个总体协方差矩阵相等的原假设

本文用到的三个excel文件如下

链接：https://pan.baidu.com/s/1oKCsLTsU99CjQciIKmPRKw 
提取码：nj03

内容来源于吴礼斌和李柏年老师主编的《MATLAB数据分析方法》，包括但不限于课本内容，记载学习记录以及自身见解，侵权联系删除。