第2章 数据描述性分析
数据描述性分析是从样本数据出发，概括分析数据的集中位置、分散程度、相互关联关系以及数据分布的正态或偏态特征等。它是进行数据分析的基础，对不同类型量纲的数据有时还要进行变换，然后再作出合理分析。本章主要介绍样本数据的基本统计量、数据的可视化、数据分布检验及数据变换等内容。
2.1 基本统计量与数据可视化
2.1.1 一维样本数据的基本统计量
描述数据的基本特征主要为集中位置和分散程度。
设从所研究的对象（即总体）X中观测得到n个观测值
x1,x2,…,xn（2.1.1）
这n个值称为样本数据,简称数据，n称为样本容量。
我们的任务就是要对样本数据进行分析，提取数据中所包含的有用信息，从而进一步对总体的特性作出推断。
1.均值、中位数、分位数与三均值
式（2.1.1）的平均值称为样本均值，记为
x=1n∑ni=1xi（2.1.2）
样本均值描述了数据取值的集中趋势（集中位置）。样本均值计算简易,但易受异常值的影响而不稳健。
将式（2.1.1）按从小到大的次序排列，排序为k的数记为x(k)(1≤k≤n)，即x(1)≤x(2)≤…≤x(n)，则
x(1),x(2),…,x(n)（2.1.3）
称为样本数据的次序统计量。
由次序统计量定义
M=xn+12,n为奇数
12xn2+xn2+1,n为偶数（2.1.4）
称M为式（2.1.1）数据的中位数。
中位数用来描述数据的中心位置的数字特征，比中位数大或小的数据个数大约为样本容量的一半。若数据的分布对称，则均值与中位数比较接近；若数据的分布为偏态，则均值与中位数差异会较大。中位数的一个显著特点是受异常值的影响较小，具有较好的稳健性。
设0≤p≤1，样本数据的p分位数定义为
Mp=x(［np］+1),np不是整数
12(x(np)+x(np+1)）,np为整数（2.1.5）
其中“\[np\]”表示np的整数部分。
显然，当p=0.5时，M0.5=M，即数据的0.5分位数等于其中位数。
一般来说，从整批数据（总体）中抽取样本数据，则整批数据中约有100p%个不超过样本数据的p分位数。在实际应用中，0.75分位数与0.25分位数比较常用，它们分别称为上、下四分位数，分别记为Q3、Q1。
一方面，虽然均值x与中位数M都是用来描述数据集中位置的数字特征的，但x是用了数据的全部信息，M只用了部分信息，因此通常情况下均值比中位数有效。另一方面，当数据有异常值时，中位数比较稳健。为了兼顾两方面的优势，人们提出三均值的概念，定义三均值如下：
M=14M0.25+12M+14M0.75（2.1.6）
按定义，三均值是上四分位数、中位数与下四分位数的加权平均，即分位数向量(M0.25,M,M0.75)与权向量（0.25,0.5,0.25）的内积。
在MATLAB中，提供了求均值、中位数、分位数等命令函数。
1）均值命令mean，其调用格式为：
m=mean(X);
其中，输入X为样本数据，输出m为样本均值。当X为矩阵时，输出为X每一列的均值向量。
2）中位数命令median，其调用格式为：
MD=median(X)；
其中，输入参数X是样本数据，输出MD为中位数。当X为矩阵时，输出为X每一列的中位数向量。
3）p分位数命令prctile，其调用格式为：
Mp=prctile(X,P)；
其中输入参数X是样本数据，P为介于0~100之间的整数，P=100*p，输出Mp为P%分位数。
注意：当样本数据X是矩阵时，上述三个命令的输出将给出X的每列数据的相对应的数值，参见例2.1.1。
4）根据分位数命令及式（2.1.6），可编写计算三均值的MATLAB程序如下：
w=\[0.25,0.5,0.25\]； %输入权向量w
SM=w*prctile(X,\[25;50;75\]);%由式(2.1.6)计算X三均值
例2.1.1 表2-1是2014年安徽省16个市森林资源情况统计数据，计算各指标均值、中位数以及三均值。
表2-1 安徽省各市森林资源情况（2014年）
地区
林业用地面积
（千公顷）
森林面积
（千公顷）
森林覆盖率
（%）
活立木总蓄积量
（万立方米）
森林蓄积量
（万立方米）
合肥市
149.05
127.33
11.13
1329.08
1002.75
淮北市
57.25
50.19
18.31
337.30
275.95
亳州市
150.93
141.65
16.62
1188.95
908.30
宿州市
279.70
255.23
25.68
1534.40
1266.15
蚌埠市
118.19
101.09
16.99
797.73
697.46
阜阳市
189.00
184.89
18.27
1225.60
1078.16
淮南市
31.37
24.27
9.39
278.74
229.39
滁州市
233.84
192.92
14.27
1608.17
1091.94
六安市
738.26
704.34
38.28
3513.36
3229.06
马鞍山市
75.09
62.32
15.39
332.15
276.64
芜湖市
134.34
102.82
17.06
533.95
438.07
宣城市
760.53
711.55
57.79
3011.54
2936.15
铜陵市
39.87
33.79
31.93
161.39
142.20
池州市
556.03
500.65
59.61
2921.23
2814.50
安庆市
619.00
575.71
37.38
3057.64
2842.54
黄山市
825.95
796.76
82.32
4742.09
4718.16
资料来源：《安徽统计年鉴2015》（http://www.ahtjj.gov.cn/tjj/tjnj/2015/2015.html）。
解：按第1章介绍的矩阵输入方法，首先将表2-1中的数据作为矩阵A输入MATLAB，然后对矩阵A调用有关计算命令，程序如下。
clear
A=\[149.05,…,4718.16\]; %输入数据，A的每一列是表2-1对应指标的样本数据
M=mean(A);%计算各指标（即各列）的均值
MD=median(A);%计算各指标的中位数
SM=\[0.25,0.5,0.25\]*prctile(A,\[25,50,75\]);%计算各指标的三均值
\[M;MD;SM\]%输出计算结果(见表2-2)
表2-2 安徽省森林资源均值、中位数与三均值（2014年）
统计量
林业用地面积
森林面积
森林覆盖率
活立木总蓄积量
森林蓄积量
均值
309.9
285.3
29.4
1660.8
1496.7
中位数
167.0
163.3
18.3
1277.3
1040.5
三均值
256.0
236.6
22.6
1489.2
1316.7
2.方差、变异系数与高阶矩
方差是描述数据取值分散性的一种度量，它是数据相对于均值的偏差平方的平均。样本数据的方差记为
s2=1n-1∑ni=1(xi-x)2=1n-1∑ni=1x2i-nx2（2.1.7）
其算术平方根称为标准差或根方差，即
s=1n-1∑ni=1x2i-nx2（2.1.8）
变异系数是描述数据相对分散性的统计量，其计算公式为
v=s/x 或者 v=s/x（2.1.9）
变异系数是一个无量纲的量，一般用百分数表示。
样本的k阶原点矩定义为
ak=1n∑ni=1xki (k=1,2,3,…)（2.1.10）
样本的k阶中心矩定义为
bk=1n∑ni=1(xi-x)k (k=1,2,3,…)（2.1.11）
在MATLAB中，计算方差命令为var，调用格式为：
S=var(X,flag);
其中，输入X为样本数据，可选项flag默认取0，输出S为方差。若flag取1，表示未修正的样本方差s20=1n∑ni=1(xi-x)2。
计算标准差命令std的调用格式为：
d=std(x, flag)；
其中，输入x是样本数据，可选项flag默认取0，输出d为标准差。若flag取1，表示未修正的样本标准差d0=s20。当输入X是矩阵时，输出X每列数据的方差或标准差。
由均值与方差命令，可编写变异系数的计算程序为：
v=std(x)./mean(x)
或者
v=std(x)./abs(mean(x))
当输入x是矩阵时，输出x每列数据的变异系数。
由均值命令mean，可编程计算k阶原点矩与中心矩，程序为：
ak=mean(x.^k) %k阶原点矩
bk=mean((x-mean(x)).^k) %k阶中心矩
MATLAB中还提供了中心矩命令moment,调用格式为：
bk=moment(X,k)；
其中，X为样本数据，k为矩的阶数。
例2.1.2 计算例2.1.1中各指标的方差、标准差和变异系数。
解：将表2-1中的数据作为矩阵A输入，然后调用有关计算命令，程序如下。
clear
A=\[149.05,…,4718.16\]; %输入原始数据（注：为节约篇幅，大部分数据用省略号表示了）
M=mean(A);%计算各指标均值
a2=mean(A.^2);%计算各指标的2阶矩
D=var(A);%计算各指标方差
SD=std(A);%计算各指标标准差
V=SD./abs(M)%计算各指标变异系数
\[D;SD;V\]%输出计算结果
将结果整理为如表2-3所示。
表2-3 安徽省森林资源的方差、标准差与变异系数（2014年）
统计量
林业用地面积
森林面积
森林覆盖率
活立木总蓄积量
森林蓄积量
方差
81232.35
74554.78
437.19
1901430.76
1877342.59
标准差
285.01
273.05
20.91
1378.92
1370.16
变异系数
0.92
0.96
0.71
0.83
0.92
3.样本的极差与四分位极差
式（2.1.1）样本数据的极大值与极小值的差称为极差，其计算公式为：
R=x(n)-x(1)（2.1.12）
极差是一种较简单的表示数据分散性的数字特征。
样本数据上、下四分位数Q3、Q1之差称为四分位极差,即
R1=Q3-Q1（2.1.13）
四分位极差也是度量数据分散性的一个重要数字特征。由于分位数对异常值有抗扰性，所以四分位极差对异常数据亦具有抗扰性。
在MATLAB中，求极差的命令为range，调用格式为：
R=range(x)；
其中，输入x是样本数据，输出R是极差。计算四分位极差的命令为iqr，调用格式为：
R1=iqr(x)；
其中，输入x是样本数据，输出R1是四分位极差。
4.异常点判别
在解决实际问题时需要对异常数据进行处理。一般判别异常值的比较简单的方法是：先计算数据的上截断点
Q3+1.5R1
与下截断点
Q1-1.5R1
其次，将数据逐个与截断点比较，小于下截断点的数据为特小值，大于上截断点的数据为特大值，两者均判为异常值。
例2.1.3 根据2013年华东各地区高校教职工数据（见表2-4），计算专任教师、行政人员、教辅人员以及工勤人员占在职教工的百分比，该百分比的极差、四分位极差以及上、下截断点。
表2-4 2013年华东各地区高校教职工数据
（单位：人）
地区
教职工数
专任教师
行政人员
教辅人员
工勤人员
上海
73361
40297
12397
8822
5693
江苏
166223
108272
22568
15056
10080
浙江
85381
56000
13712
7878
3430
安徽
76178
54903
8105
5947
4618
福建
64744
42905
9719
5909
3314
江西
74396
52434
8739
5693
4054
山东
142240
98685
17617
11989
8639
数据来源：《中国统计年鉴2014》（http://www.stats.gov.cn/tjsj/ndsj/2014/zk/html/Z2131C.xls）。
解：将表2-4中的数据作为矩阵A输入，然后调用有关计算命令，程序如下：
clear
A=\[73361 40297 12397 8822 5693
166223 108272 22568 15056 10080
85381560001371278783430
7617854903810559474618
6474442905971959093314
7439652434873956934054
1422409868517617119898639\]; %输入表2-4中数据
B=A(:,2:5)./\[A(:,1)*ones(1,4)\];%计算百分比
R=range(B);%计算百分比极差
R1=iqr(B);%计算四分位极差
XJ=prctile(B,\[25\])-1.5*R1;%计算下截断点
SJ=prctile(B,\[75\])+1.5*R1;%计算上截断点
ycz1=B<ones(7,1)*XJ;%小于下截断点的数据
ycz2=B>ones(7,1)*SJ;%大于上截断点的数据
由程序运行结果可以得知：上海市专任教师占在职教工的百分比小于其他省份，教辅人员以及工勤人员占在职教工的百分比大于其他省份，可认为上海市这3项指标为异常值。
5.偏度与峰度
偏度与峰度是样本数据分布特征和正态分布特征比较而引入的概念。偏度是用于衡量分布的非对称程度或偏斜程度的数字特征。样本数据的偏度定义为：
sk=b3b3/22=1n∑ni=1(xi-x)31n∑ni=1(xi-x)23（2.1.14a）
或样本数据修正的偏度定义为：
sk=n2b3(n-1)(n-2)s3（2.1.14b）
其中b2、b3、s分别表示样本的2、3阶中心矩与标准差。
当sk>0时，称数据分布右偏，此时均值右边的数据比均值左边的数据更散，分布的形状是右长尾的；当sk<0时，称数据分布左偏，均值左边的数据比均值右边的数据更散，分布的形状是左长尾的；当sk接近于0时，称分布无偏倚即认为分布是对称的。正态分布的样本数据的偏度接近于0，当样本数据的偏度与零相差较大，则可初步拒绝样本数据来自于正态分布总体。
一般有：数据分布右偏时，算术平均数>中位数>众数；左偏时相反，众数>中位数>算术平均数。
峰度是用来衡量数据尾部分散性的指标。样本数据的峰度定义为：
ku=b4b22-3=1n∑ni=1(xi-x)41n∑ni=1(xi-x)22-3（2.1.15a）
或样本数据修正的峰度定义为：
ku=n2(n+1)b4(n-1)(n-2)(n-3)s4-3(n-1)2(n-2)(n-3)（2.1.15b）
其中b4、s分别表示样本数据的4阶中心矩与标准差。
当数据的总体分布是正态分布时，峰度近似为0；与正态分布相比较，当峰度大于0时，数据中含有较多远离均值的极端数值，称数据分布具有平峰厚尾性；当峰度小于0时，表示均值两侧的极端数值较少，称数据分布具有尖峰细尾性。在金融时间序列分析中，通常需要研究数据是否为尖峰、厚尾等特性。
在MATLAB中，计算样本数据的偏度命令为skewness，调用格式为：
sk=skewness(x,flg)；
其中输入x是样本数据，flg取0或1。当flg取0时，是修正的偏度；当flg取1时，按式（2.1.14b）计算偏度，系统默认flg=1。当x是矩阵时，输出sk为数组，其第i个元素是x的第i列数据的偏度。
与样本数据的峰度有关的命令为kurtosis，调用该命令计算峰度的程序为：
ku=kurtosis (x, flg)-3；
其中输入x是样本数据，flg取0或1。当flg取0时，是修正的峰度，当flg取1时，按式（2.1.15b）计算峰度，系统默认flg=1。输出ku为峰度，当x是矩阵时，ku为数组，其第i个元素是x的第i列数据的峰度。（注意:命令kurtosis给出的结果与我们峰度的定义相差常数3。）
例2.1.4 “中国银行”股票的交易数据（2015年9月1日到2016年3月18日）（在第1章中收集保存在文件'sh601988.xls'中），计算股票收盘价、最高价、最低价、开盘价，以及涨跌幅数据的偏度、峰度。
解：在MATLAB命令窗口中输入程序如下。
clear
zgyh=xlsread('sh601988.xls','sheet1') %读取数据文件sh601988.xls
pd=skewness(zgyh,0);%计算zgyh每列数据的偏度
fd=kurtosis(zgyh,0)-3;%计算zgyh每列数据的峰度
\[pd;fd\]%输出计算结果
subplot(2,3,1),histfit(zgyh(:,1)),title('收盘价')%作收盘价直方图
subplot(2,3,2),histfit(zgyh(:,2)),title('最高价')%作最高价直方图
subplot(2,3,3),histfit(zgyh(:,3)),title('最低价')%作最低价直方图
subplot(2,3,4),histfit(zgyh(:,4)),title('开盘价')%作开盘价直方图
subplot(2,3,5),histfit(zgyh(:,5)),title('涨跌幅')%作涨跌幅直方图
从表2-5中可得：2015年9月1日至2016年3月18日“中国银行”股票收盘价、最高价、最低价、开盘价以及涨跌幅数据偏度均小于0，因此数据分布均左偏；涨跌幅峰度大于0，且较大，数据具有平峰厚尾性，其余的峰度小于0，具有尖峰细尾性。总之可认为数据总体不服从正态分布。从各个指标数据的直方图（如图2-1所示）也可直观看到分布呈左偏态。
表2-5 “中国银行”股票收盘价、最高(低)价、开盘价、涨跌幅偏度与峰度
统计量
收盘价
最高价
最低价
开盘价
涨跌幅
偏度
-0.5784
-0.5386
-0.6251
-0.6352
-0.0874
峰度
-1.1048
-1.0363
-1.1331
-1.0840
2.1097
图2-1 “中国银行”股票各指标直方图
2.1.2 多维样本数据的统计量
1.多维样本数据的统计量
多维样本数据的统计量主要有：样本均值向量、样本协方差矩阵、样本相关系数矩阵等。
设总体为p维向量G=(X1,X2,…,Xp)，从中抽取样本容量为n的样本，第i个样本观测值为(xi1,xi2,…,xip)(i=1,2,…,n)，记
X=x11x12…x1p
x21x22…x2p
xn1xn2…xnp（2.1.16）
称X为样本数据矩阵。为了方便起见，X的第j个列向量记为Xj=(x1j,x2j,…,xnj)T。
显然，X的第j个列向量是Xj的n个观测数据。通常由样本数据矩阵X出发，构造下列统计量来分析总体的特征。
1）样本均值向量。记Xj的观测值（即X中的第j列）的均值为
xj=1n∑ni=1xij (j=1,2,…,p)（2.1.17）
称x=(x1,x2,…,xp)T为p维样本均值向量。
2）样本协方差矩阵。记
sjk=1n-1∑ni=1(xij-xj)(xik-xk) (j,k=1,2,…,p)（2.1.18）
称sjk为Xj与Xk的样本协方差，或称为样本数据矩阵X的第j列与第k列的协方差。记
S=s11s12…s1p
s21s22…s2p
sp1sp2…spp（2.1.19）
称S为样本协方差矩阵。
显然，Xj的方差sjj，即
sjj=1n-1∑ni=1(xij-xj)2 (j=1,2,…,p)（2.1.20）
3)样本相关系数矩阵。X的第j列与第k列的相关系数记为
rjk=sjksjjskk (j,k=1,2,…,p)
又记
R=r11r12…r1p
r21r22…r2p
rp1rp2…rpp=1r12…r1p
r211…r2p
rp1rp2…1（2.1.21）
称R为样本相关系数矩阵。
不难验证，样本相关系数矩阵与协方差矩阵存在如下关系：
R=DTSD（2.1.22）
其中D=diag(s-1/211,s-1/222,…,s-1/2pp)。
4）样本标准化矩阵。令
x*ij=xij-xjsjj(i=1,2,…,n;j=1,2,…,p)（2.1.23）
记X*=(x*ij)n×p（2.1.24）
称为样本矩阵X的标准化矩阵。
可以证明：X*的协方差矩阵S*等于X的相关系数矩阵R，即S*=R。
5）Rc矩阵。
矩阵X的第j列与第k列的Rc系数定义为：
rcjk=2［Xj,Xk］［Xj,Xj］+［Xk,Xk］（2.1.25）
其中［Xj,Xk］=∑ni=1xijxik(j,k=1,2,…,p)，称矩阵(rcjk)p×p为矩阵X的Rc矩阵，记为Rc(X)，即
Rc(X)=(rcjk)p×p
由定义式（2.1.25），显然rcjj=1(j=1,2,…,p)，rcjk≤1(j,k=1,2,…,p)。
可以证明：对于矩阵X，有Rc(X*)=R，即X的标准化矩阵的Rc矩阵等于其相关系数矩阵。
协方差矩阵与量纲有关，相关系数矩阵R以及Rc矩阵与量纲无关，这一点在今后的判别分析中值得注意。
协方差矩阵S、相关系数矩阵R与Rc矩阵都是实对称非负定矩阵。
2.多维样本数据的MATLAB命令
在MATLAB中，计算样本协方差矩阵命令为cov，调用格式为：
S=cov(X)；
当X为向量时,S表示X的方差；当X为矩阵时,S为X的协方差矩阵，即S的对角线元素是X每列的方差,S的第i行第j列元素为X的第i列和第j列的协方差值。
计算样本相关系数矩阵命令为corr，调用格式为：
R= corr (X)；
其中,X为样本矩阵，输出R的对角线元素是1，R的第i行第j列元素为X的第i列和第j列的相关系数。
计算X的标准化矩阵命令为zscore,调用格式为：
Z= zscore (X)；
其中X为样本矩阵，输出Z是标准化矩阵。
MATLAB中没有计算Rc矩阵的命令，因此根据Rc矩阵的定义，可编写计算Rc矩阵的程序如下。
%输入样本数据矩阵X
X=\[data\];
%计算RC矩阵
for i=1: size(X,2)
for j=1: size(X,2)
RC(i,j)=2*dot(X(:,i),X(:,j))./\[sum(X(:,i).^2)+ sum(X(:,j).^2)\];
end
end
RC%输出RC（X）
例2.1.5 中国银行（股票代码：601988）、交通银行（601328）、工商银行（601398）、建设银行（601939）与农业银行（601288）股票2016年1月4日至2016年3月31日的收盘价见表2-6（数据保存于文件“yhgspj.xls”中），求5支股票收盘价间的协方差、相关系数、Rc系数矩阵以及收盘价的标准化矩阵。
表2-6 5支银行股的收盘价
（单位：元）
日期
中国银行
(601988)
交通银行
(601328)
工商银行
(601398)
建设银行
(601939)
农业银行
(601288)
2016/1/4
3.87
6.07
4.45
5.6
3.12
2016/1/5
3.86
6.16
4.47
5.58
3.16
2016/1/6
3.89
6.21
4.51
5.64
3.18
2016/1/7
3.79
5.94
4.43
5.49
3.11
…
…
…
…
…
…
2016/3/30
3.4
5.58
4.32
4.9
3.2
2016/3/31
3.4
5.57
4.29
4.85
3.2
数据来源：上海证券交易所网站（http://www.sse.com.cn/）。
解：
clear
a=xlsread('yhgspj.xls');%读取收盘价数据
a1=flipud(a);%矩阵上下翻转，按交易日期顺序
S=cov(a1);%计算协方差矩阵
R=corr(a1);%计算相关系数矩阵
Z=zscore(a1);%数据标准化
for i=1: size(a1,2)%计算Rc系数矩阵
for j=1: size(a1,2)
RC(i,j)=2*dot(a1(:,i),a1(:,j))./\[sum(a1(:,i).^2)+ sum(a1(:,j).^2)\];
end
end
RC
输出结果：
S =
0.0364 0.0484 0.0242 0.0485 0.0127
0.0484 0.0681 0.0342 0.0626 0.0199
0.0242 0.0342 0.0208 0.0284 0.0134
0.0485 0.0626 0.0284 0.0710 0.0119
0.0127 0.0199 0.0134 0.0119 0.0110
R =
1.0000 0.9715 0.8813 0.9542 0.6362
0.9715 1.0000 0.9106 0.9009 0.7250
0.8813 0.9106 1.0000 0.7405 0.8846
0.9542 0.9009 0.7405 1.0000 0.4236
0.6362 0.7250 0.8846 0.4236 1.0000
RC =
1.0000 0.8929 0.9763 0.9348 0.9946
0.8929 1.0000 0.9651 0.9933 0.8525
0.9763 0.9651 1.0000 0.9877 0.9531
0.9348 0.9933 0.9877 1.0000 0.8996
0.9946 0.8525 0.9531 0.8996 1.0000
从相关系数矩阵可知：农业银行与建设银行在此期间相关系数最小，而中国银行与交通银行相关系数最大。
2.1.3 样本数据可视化
1.可视化
数据可视化是指数据的图形表示。借助几何图形可形象地说明数据的特征与分布情况。常用的图形有条形图、直方图、盒图、阶梯图和火柴棒图等。
1）条形图。条形图是用宽度相同的直线条的高低或长短来表示统计指标数值的大小。条形图根据表现资料的内容可分为单式条形图、复式条形图和结构条形图。单式条形图反映统计对象随某一因素变化而改变的情况。复式条形图可以反映统计对象随两个因素变动而变动的情况。结构条形图则反映不同统计对象内部结构的变化情况。
在MATLAB中，绘制条形图命令bar的调用格式为：
① bar(X)；
② bar(x,Y)；
①作样本数据X的条形图；②x的元素在横坐标轴上按从小到大排列，作Y和x对应的条形图。
2）直方图。将观测数据的取值范围分为若干个区间,计算落在每个区间的频数或频率。在每个区间上画一个矩形，以估计总体的概率密度。
在MATLAB中，绘制直方图命令hist的调用格式为：
① hist(x,n)； %作数据x的直方图，其中n表示分组的个数，缺省时n=10
② \[h,stats\] = cdfplot(x)；%作数据x的经验分布函数图，stats给出数据的最大值、最小值、中位数、平均值和标准差。
附加有正态密度曲线的直方图命令histfit的调用格式为：
① histfit(X)%X为样本数据向量，返回直方图和正态曲线
② histfit(X,nbins)%nbins指定bar的个数，缺省时为X中数据个数的平方根
3）盒图。盒图由五个数值点组成：最小值、下四分位数、中位数、上四分位数、最大值。中间的盒子从Q1延伸到Q3，盒子里的直线标示出中位数的位置，盒子两端有直线往外延伸到最小数与最大数。
在MATLAB中，绘制盒图命令boxplot的调用格式为：
boxplot(X)；%产生矩阵X的每一列的盒图和“须”图，“须”是从盒的尾部延伸出来，并表示盒外数据长度的线，如果“须”的外面没有数据，则在“须”的底部有一个点
4）阶梯图。命令stairs的调用格式为：
stairs(x)；%作数据x的阶梯图
5）火柴棒图。命令stem的调用格式为：
stem(x)%作离散数据序列x的火柴棒图
例2.1.6 随机生成150个服从标准正态分布随机数，将这些数据作为样本数据，分别作出样本数据的柱形图、直方图、阶梯图、火柴棒图等图形。
解：编写程序如下。
clear
x = random('normal',0,1,\[1,150\]); %产生服从标准正态分布随机数150个
bar(x)%作柱形图（如图2-2所示）
figure(2);hist(x)%作直方图（如图2-3所示）
figure(3);stairs(x)%作阶梯图（如图2-4所示）
figure(4);stem(x)%作火柴棒图（如图2-5所示）
图2-2 柱形图
图2-3 直方图
图2-4 阶梯图
图2-5 火柴棒图
读者还可作出盒图与附加有正态密度曲线的直方图等。
2.二维与三维数据可视化
1）绘制散点图的命令scatter与scatter3的调用格式为：
scatter(x,y)；
其中x是横坐标向量，y是纵坐标向量，输出平面散点图。
scatter3(x,y,z)；
其中x、y、z分别是横、纵、竖坐标向量，输出空间散点图。
2）绘制曲面图的命令mesh与surf的调用格式为：
mesh(X,Y,Z)
或者
surf(X,Y,Z)
其中Z是对应(X,Y)处的函数值Z=f(X,Y)，\[X,Y\]是由命令meshgrid生成的数据点矩阵，即\[X,Y\]=meshgrid(x,y)，输入向量x为Oxy平面上矩形定义域的矩形分割线在x轴上的坐标，向量y为Oxy平面上矩形定义域的矩形分割线在y轴上的坐标。矩阵X为Oxy平面上矩形定义域的矩形分割点的横坐标值矩阵，X的每一行是向量x，且X的行数等于y的维数；矩阵Y为Oxy平面上矩形定义域的矩形分割点的纵坐标值矩阵，Y的每一列是向量y，且Y的列数等于x的维数。
例2.1.7 设总体(X,Y)服从二维正态分布N(2,1;3,3;3/2)，生成100对服从该分布随机数据对(xi,yi)，将这些数据作为样本数据，作出样本数据的散点图。再根据二维正态分布的密度函数，绘制密度曲面图形。
解：随机生成服从二维正态分布的数据的命令mvnrnd，其调用格式为：
X=mvnrnd(mu,sigma,n)；
其中mu是均值向量，sigma是协方差矩阵，n是数据个数，输出X是与协方差矩阵同阶的随机数据矩阵。
已知二维正态分布中的参数μ1=2,σ21=1;μ2=3,σ22=3;ρ=3/2，所以均值向量为μ=(2,3),协方差矩阵为Σ=σ21ρσ1σ2
ρσ1σ2σ22=11.5
1.53，编写程序如下：
clear
mu = \[2 3\]; %输入均值向量
sa = \[1 1.5; 1.5 3\];%输入协方差矩阵
Nr = mvnrnd(mu,sa,100);%随机生成n=100的样本数据
scatter(Nr(:,1),Nr(:,2),'*');%作样本数据平面散点图，如图2-6所示
%绘制密度曲面
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)%绘制密度曲面图形，如图2-7所示
输出图形结果如下：
图2-6 样本数据的散点图
图2-7 服从N(2,1;3,3;3/2)分布的密度曲面图
由图形2-6可以看出，散点图位于平面上的一个椭圆状区域内，不同的相关系数对应的椭圆状区域形状不同，相关系数越接近于1，椭圆越扁长，可以利用这一图形特征初步说明数据是否来自正态总体。
3.正态概率图与Q-Q图
1）正态概率图。正态概率图用于正态分布的检验，其横坐标是样本数据的分位数，纵坐标是标准正态分布的α分位数对应的概率α。设总体X服从正态分布N(μ,σ2)，来自总体的样本为x1,x2,…,xn，其次序统计量x(1)≤x(2)≤…≤x(n)，则X的经验分布函数为Fi=P{X≤x(i)}，Fi的估计取为Fi=i-0.375n+0.25。在平面上的n个点
(x(i),Φ-1(Fi))(i=1,2,…,n)
的散点图称为正态概率图。此图形的纵坐标不以Φ-1=Φ-1(F)为刻度，而标以相应的F值，这种坐标系也称为概率纸坐标系（如图2-8所示）。可以证明：若样本是来自正态总体，则正态概率图呈现一条直线形状。
在MATLAB中，绘制正态概率图的命令为normplot，调用格式为：
normplot(X)；
其中X为样本数据。对于输出的正态概率图，每一个样本数据对应图上的一个“+”号，图中有一条红色参考直线，若图中的“+”号都集中在这条参考线附近，说明样本数据近似服从正态分布，否则偏离参考线的“+”号越多，说明样本数据不服从正态分布。
同理，总体是非正态分布时，也可以类似地绘制概率图。如MATLAB系统中，绘制威布尔分布概率图的命令为weibplot，调用格式为：
weibplot(X)；
其中，当输入X为向量时，显示威布尔（Weibull）分布概率图。如果样本数据点基本散布在一条直线上，则表明数据服从该分布，否则拒绝该分布。
例2.1.8 对于例2.1.7模拟的样本数据Nr，分别作出两个分量的正态概率图，从图形直观检验各分量是否服从正态分布。
解：编写程序如下（接例2.1.6程序）。
figure(3)
subplot(1,2,1),normplot(Nr(:,1))%分量x正态概率图
subplot(1,2,2),normplot(Nr(:,2))%分量y的正态概率图
从图2-8可以看出，两个分量的正态概率图呈现一条直线形状，所以样本中的每个分量都服从正态分布。事实上，数学理论上已证明：二元正态分布的边缘分布仍为正态分布。正态概率图反映的结果与理论相一致。
2）Q-Q图。Q-Q图可以用于检验样本数据是否服从指定的分布，是样本数据的分位数与所指定分布的分位数之间的关系曲线图。通常情况下，一个坐标轴表示样本分位数，另一个坐标轴表示指定分布的分位数。每一个样本数据对应图上的一个“+”号，图中有一条参考直线，若图中的“+”号都集中在这条参考线附近，说明样本数据近似服从指定的分布；否则偏离参考线的“+”号越多，说明样本数据越不服从指定的分布。
例如：设总体服从正态分布N(μ,σ2)，来自总体的样本为x1,x2,…,xn，其次序统计量x(1)≤x(2)≤…≤x(n)，则平面上n个点
Φ-1i-0.375n+0.25,x(i)(i=1,2,…,n)（2.1.26）
的散点图称为Q-Q图，其中Φ-1()为标准正态分布函数的反函数。
图2-8 x分量（左）与y分量（右）的正态概率图
可以证明，若样本确是来自正态总体,则散点在直线y=σx+μ附近，即Q-Q图大致呈现一条直线形状；当样本来自其他分布总体时，样本Q-Q图将是弯曲的。这样，利用Q-Q图可以直观地作正态性检验，即若Q-Q图近似一条直线时，则可认为样本数据来自正态总体。
对于其他类型的分布，也有相应的Q-Q图，其中散点的横坐标为该分布的对应分位数，以此可判断数据是否近似服从该类型的分布。
在MATLAB中，绘制Q-Q图的命令为qqplot,调用格式为：
qqplot(X)；
或
qqplot(X,PD)；
其中X为样本数据，PD是由fitdist命令指定的分布类，其指定调用方式为：
PD=fitdist(X,distname)；
这里的distname是指定的分布类名称。常用的distname有'Beta'（贝塔分布）、'Binomial'（二项分布）、'Exponential'（指数分布）、'Normal'（正态分布）、'Weibull'（威布尔分布）等。qqplot中省略PD时默认为标准正态分布。
例2.1.9 对于例2.1.7模拟的样本数据Nr，分别作出两个分量的Q-Q图，从图形直观检验各分量是否服从正态分布。
解：编写程序如下（接例2.1.6程序）。
figure(4)
subplot(1,2,1),qqplot(Nr(:,1)),grid on%分量x正态Q-Q图
subplot(1,2,2),qqplot(Nr(:,2)),grid on%分量y的正态Q-Q图
从图2-9可以看出，两个分量的Q-Q图呈现一条直线形状，所以样本中的每个分量都服从正态分布。
图2-9 x分量（左）与y分量（右）的Q-Q图
例2.1.10 模拟生成服从自由度为8的卡方分布的样本数据300个，记为c1，分别绘制c1的正态概率图与卡方分布的Q-Q图，从图形直观检验该数据是否服从指定的分布。
解：程序如下。
clear
s=rng; rng(s); %保持生成样本不变，即程序每次运行时生成的随机数相同
c1=chi2rnd(8,\[300,1\]); c2=sort(c1); %模拟生成卡方分布样本
plot(c2,chi2pdf(c2,8), '+-');%绘制卡方分布的密度曲线（如图2-10所示）
title('卡方分布的密度曲线') ;legend('自由度n=8');
grid on
figure
pd=makedist('Gamma','a',4,'b',0.5)%创建参数a=4,b=0.5的伽马分布
subplot(1,2,1),normplot(c1);%绘制样本的正态概率图（如图2-11 a所示）
subplot(1,2,2),qqplot(c1,pd);%按指定分布绘制样本的Q-Q图（如图2-11 b所示）
grid on
图2-10 样本的卡方分布的密度曲线
图2-11 样本的正态概率图与Q-Q图
从图2-10可看出样本数据分布是偏态的。图2-11a即正态概率图两端偏离参考直线明显，表明数据不服从正态分布。图2-11b即Q-Q图呈现在参考直线附近，可初步判定数据服从指定的伽马分布，程序中指定的伽马分布参数a=4、b=0.5，此即是自由度为8的卡方分布（注：自由度为n的卡方分布是参数a=n/2,b=0.5的伽马分布）。
