被除名的飓风:matlab 实验05数据的统计分析(终于找到了)

实验五  数据的统计分析
一、问题背景与实验目的
二、相关函数（命令）及简介
三、实验内容
四、自己动手
五、附录
一、问题背景与实验目的
在日常生活中我们会在很多事件中收集到一些数据（比如：考试分数、窗口排队人数、月用电量、灯泡寿命、测量误差、产品质量、月降雨量等数据），这些数据的产生一般都是随机的．这些随机数据乍看起来并没有什么规律，但通过数理统计的研究发现：这些随机数还是符合着某种分布规律的，这种规律被称为统计规律．
本实验旨在通过对概率密度函数曲线的直观认识、对数据分布的形态猜测、对某些概率分布的密度函数的参数估计（以正态为例）以及进行简单的正态假设检验，来揭示生活中的随机数据的一些统计规律．
二、相关函数（命令）及简介
1. 概率密度函数pdf系列．以normpdf( )为例，调用格式：
y=normpdf(x, mu,sigma)，
计算参数为mu和sigma的样本数据x的正态概率密度函数．参数sigma必须为正．其中：mu为均值，sigma为标准差．
2. 参数估计fit系列．以normfit( )为例，调用格式：
[muhat, sigmahat, muci, sigmaci] = normfit(x, alpha)，
对样本数据x进行参数估计，并计算置信度为100（1－alpha）%的置信区间．如alpha=0.01时，则给出置信度为99％的置信区间．不写明alpha，即表示alpha取0.05．
3．load( )函数．调用格式：
S = load('数据文件')
将纯数据文件（文本文件）中的数据导入Matlab，S 是双精度的数组，其行数、列数与数据文件相一致．
4. hist(x, m)函数：画样本数据x的直方图，m为直方图的条数，缺省值为10．
5. tabulate( )函数：绘制频数表．返回table矩阵，第一列包含x的值，第二列包含该值出现次数，最后一列包含每个值的百分比．
6．ttest(x,m,alpha) 函数：假设检验函数．此函数对样本数据x进行显著性水平为alpha的t假设检验，以检验正态分布样本x（标准差未知）的均值是否为m．h=1表示拒绝零假设，h=0表示不能拒绝零假设．
7．normplot(x)或weibplot(x) 函数：统计绘图函数，进行正态分布检验．
研究表明：如果数据是来自一个正态分布，则该线为一直线形态；如果它是来自其他分布，则为曲线形态．
完全类似地可探索以下一系列函数的用法与作用：
8．累积分布函数cdf系列，如：normcdf( )．
9．逆累积分布函数inv系列，如：norminv( )．
10．随机数发生函数rnd系列，如：normrnd( )．
11．均值与方差函数stat系列，如：normstat( )．
三、实验内容
1. 常见的概率分布的密度函数及其图形
1）常见概率分布的密度函数（20个，打√的10个将在后面作介绍）
序号
中文函数名
英文函数名
英文简写
备注
1
Beta分布
Beta
beta
2
二项分布
Binomial
bino
√
3
卡方分布
Chisquare
chi2
√抽样
4
指数分布
Exponential
exp
√
5
F分布
F
f
√抽样
6
Gamma分布
Gamma
gam
7
几何分布
Geometric
geo
√
8
超几何分布
Hypergeometric
hyge
9
对数正态分布
Lognormal
logn
10
负二项式分布
Negative Binomial
nbin
11
非中心F分布
Noncentral F
ncf
12
非中心t分布
Noncentral t
nct
13
非中心卡方分布
Noncentral Chi-square
ncx2
14
正态分布
Normal
norm
√
15
泊松分布
Poisson
poiss
√
16
瑞利分布
Rayleigh
rayl
17
T分布
T
t
√抽样
18
均匀分布
Uniform
unif
√
19
离散均匀分布
Discrete Uniform
unid
√
20
Weibull分布
Weibull
weib
2）常见概率分布的密度函数文字说明与图形演示：
（A）常见连续分布的密度函数
（1）正态分布
若连续型随机变量的密度函数为：

则称为服从正态分布的随机变量，记作．特别地，称时的正态分布为标准正态分布，其概率分布的密度函数参见图1．一个非标准正态分布的密度函数参见图2中的虚线部分（）．
正态分布是概率论与数理统计中最重要的一个分布，高斯(Gauss)在研究误差理论时首先用正态分布来刻画误差的分布，所以正态分布又称高斯分布．一个变量如果是由大量微小的、独立的随机因素的叠加效果，那么这个变量一定是正态变量．比如测量误差、产品质量、月降雨量等都可用正态分布描述．
x=-8:0.1:8;
y=normpdf(x, 0, 1);
y1=normpdf(x, 1, 2);
plot(x, y, x, y1, ':' );
 
图1 标准正态分布             图2 标准正态与非标准正态
（2）均匀分布（连续）
若随机变量的密度函数为

则称服从区间上的均匀分布（连续），记作，其概率分布的密度函数见参见图3．
均匀分布在实际中经常使用，譬如一个半径为的汽车轮胎，因为轮胎上的任一点接触地面的可能性是相同的，所以轮胎圆周接触地面的位置是服从上的均匀分布，这只要看一看报废轮胎四周磨损程度几乎是相同的就可明白均匀分布的含义了．
x=-10:0.01:10;r=1;
y=unifpdf(x, 0, 2*pi*r);
plot(x, y);

图3均匀分布（连续）            图4　指数分布
（3）指数分布
若连续型随机变量的密度函数为：
  其中，
则称为服从参数为的指数分布的随机变量，记作．
在实际应用问题中，等待某特定事物发生所需要的时间往往服从指数分布．如某些元件的寿命；某人打一个电话持续的时间；随机服务系统中的服务时间；动物的寿命等都常假定服从指数分布．
指数分布的重要性还在于它是具有无记忆性的连续型随机变量．即：设随机变量服从参数为的指数分布，则对任意的实数，有

其概率分布的密度函数参见见图4．
x=0:0.1:30;
y=exppdf(x, 4);
plot(x, y)
（B）常见离散分布的密度函数
（4）几何分布
在一个贝努里实验中，每次试验成功的概率为，失败的概率为，设试验进行到第次才出现成功，则的分布列为：

容易看到是几何级数的一般项，于是人们称它为几何分布，其概率分布的密度函数参见图5．
x=0:30;
y=geopdf(x, 0.5);
plot(x, y)

图5 几何分布                  图6 二项分布
（5）二项分布
如果随机变量的分布列为：

则这个分布称为二项分布，记为．当时的二项分布又称为0-1分布，分布律为

0
1



一般的二项分布的密度函数参见图6.
x=0:50;
y=binopdf(x, 500, 0.05);
plot(x, y);
（6）泊松(Poisson)分布
泊松分布是1837年由法国数学家泊松（Poisson S.D.1781-1840）首次提出的，其概率分布列是：

记为，其概率分布的密度函数参见图7．
泊松分布是一种常用的离散分布，它与单位时间（或单位面积、单位产品等）上的计数过程相联系，譬如：单位时间内，电话总机接到用户呼唤次数；1平方米内，玻璃上的气泡数；一铸件上的砂眼数；在单位时间内，某种放射性物质分裂到某区域的质点数等等．
x=0:50;
y=poisspdf(x, 25);
plot(x, y);
注：对比二项分布的概率密度函数图可以发现，当二项分布的与泊松分布充分接近时，两图拟合程度非常高（图6与图7中的），直观地验证了泊松定理（泊松分布是二项分布的极限分布），请对比图6与图7．

图7 泊松分布             图8 均匀分布（离散）
（7）均匀分布（离散）
如果随机变量的分布列为：

则这个分布称为离散均匀分布，记为，其概率分布的密度函数参见图8．
n=20;
x=1:n;
y=unidpdf(x, n);
plot(x, y, 'o-' );
（C）三大抽样分布的密度函数
（8）分布
设随机变量相互独立，且同服从正态分布，则称随机变量服从自由度为的分布，记作，亦称随机变量为变量．其概率分布的密度函数参见图9、图10，分布的密度函数解析式参见本章的附录表格．
x=0:0.1:20;                                                   x=0:0.1:20;
y=chi2pdf(x, 4);                                                  y=chi2pdf(x, 10);
plot(x, y);                                                      plot(x, y)

图9分布             图10分布
（9）分布
设随机变量，，且与相互独立，则称随机变量

服从自由度为的分布，记作，其概率分布的密度函数参见图11，即，分布的密度函数解析式参见本章的附录表格．
x=0.01:0.1:8.01;
y=fpdf(x, 4, 10);
plot(x, y)

图11分布             图12 分布
（10）分布
设随机变量，且与相互独立，则称随机变量

服从于自由度为的分布，记作，其概率分布的密度函数参见图12，即．分布的密度函数解析式参见本章的附录表格．
细心的读者可能已经发现，图12的分布图与图1、图2的正态分布十分相似．可以证明：当时，分布趋于标准正态分布．
x=-6:0.01:6;
y=tpdf(x, 4);
plot(x, y)
2．对给定数据画频数直方图（Histogram）或频数表（Frequency Table）
假定有若干个给定的数据集，它们满足上述10种分布之一，我们现在的任务就是利用画频数直方图等手段，确定它们到底服从哪一类分布．
例1：某一次书面考试的分数罗列如下，试画频数直方图．
鉴于数据的数量较大（包含有120个数据），可以先在一个文本文件中输入，保存为data1.txt．
75    69   100    80    70    74    78    59    72    73
63    79    69    81    62    87    80    66    86    75
70    85    85    64    78    65    69    67    78    72
60    50    57    83    77    79    78    74    67    83
71    67    71    74    84    74    83    75    73    74
60    91    65    69    80    63    86    67    73    80
74    68    72    80    95    61    77    85    82    71
80    76    83    69    87    76    72    69    66    86
74    87    59    81    88    75    83    71    77    81
88    67    67    76    71    76    79    79    90    62
80    85    81    75    72    57    94    91    83    78
66    74    79    74    82    79    87    76    81    68
x=load('data1.txt');
x=x(:);
hist(x)
结果参见图13．从图形形态上来看，图13较为接近图2所示的正态分布．

图13  例1的频数直方图             图14  例2的频数直方图
例2：某一次上机考试的分数罗列如下(data2.txt，包含有130个数据)，试画频数直方图．
51    70    95    91    70    83    83    96    66    61
79    79    57    85    95    83    63    71    71    72
91    60    69   100    67    87    72    50    60    63
87    98    71    74    96    55    83    67    92    78
56    62    77    79    84    55    59    61    93    56
82    61    88    97    98    95    73    79    81    87
56    92    53    57    93    89    77    89    56    92
99    86    68    57    91    57    81    65    80    99
79    95    79    86    74    56    70    61    72    81
57    75    98    89    69    61    71    77    72    78
70    73    67    59    62    86    84    93    82    80
90    94    84    89    80    67    97    73    80    94
69    64    51    51    92    62    52    86    67    97
x=load('data2.txt');
x=x(:);
hist(x)
结果参见图14．图14看上去很接近图8所示的均匀分布（离散）．
例3：以下给出上海1998年来的月降雨量的数据(data3.txt，包含有98个数据)：
1184.4  1113.4  1203.9  1170.7  975.4   1462.3  947.8
1416.0  709.2   1147.5  935    1016.3  1031.6  1105.7
849.9   1233.4  1008.6  1063.8  1004.9  1086.2  1022.5
1330.9  1439.4  1236.5  1088.1  1288.7  1115.8  1217.5
1320.7  1078.1  1203.4  1480.0  1269.9  1049.2  1318.4
1192.0  1016.0  1508.2  1159.6  1021.3  986.1   794.7
1318.3  1171.2  1161.7  791.2   1143.8  1602.0  951.4
1003.2  840.4   1061.4  958.0   1025.2  1265.0  1196.5
1120.7  1659.3  942.7   1123.3  910.2   1398.5  1208.6
1305.5  1242.3  1572.3  1416.9  1256.1  1285.9  984.8
1390.3  1062.2  1287.3  1477.0  1011.9  1217.7  1197.1
1143.0  1018.8  1243.7   909.3  1030.3  1124.4  811.4
820.9   1184.1  1107.5   991.4   901.7  1176.5  1113.5
1272.9  1200.3  1508.7   772.3   813.0  1392.3  1006.2
x=load('data3.txt');
x=x(:);
hist(x)
结果参见图15．图15看上去很接近图10所示的分布．

图15  例3的频数直方图             图16  例4的频数直方图
在重复数据较多的情况下，我们也可以利用Matlab自带的函数tabulate( )产生频数表，并以频数表的形式来发掘数据分布的规律．
例4：给出以下数据：(data4.txt，含有46个数据)
2  3  6  4  1  5  1  2  3  1  4  2  3  1  3  3  2  3  1  6  4  6  4
6  5  4  3  6  4  3  3  3  3  4  4  5  6  2  1  2  3  4  5  6  5  4
则：
x=load('data4.txt');
x=x(:);
tabulate(x)
hist(x, 6)
Value    Count    Percent
1        6     13.04%
2        6     13.04%
3       12     26.09%
4       10     21.74%
5        5     10.87%
6        7     15.22%
结果参见图16．图16看上去好象没有什么规律可循．
例5：现累积有100次刀具故障记录，当故障出现时该批刀具完成的零件数如下：(data5.txt)
459  362  624  542  509  584  433  748   815  505
612  452  434  982  640  742  565  706   593  680
926  653  164  487  734  608  428  1153  593  844
527  552  513  781  474  388  824  538   862  659
775  859  755   49  697  515  628  954   771  609
402  960  885  610  292  837  473  677   358  638
699  634  555  570   84  416  606  1062  484  120
447  654  564  339  280  246  687  539   790  581
621  724  531  512  577  496  468  799   544  645
764  558  378  765  666  763  217  715   310  851
x=load('data5.txt');
x=x(:);
hist(x)       %%结果参见图17，很象图2所示的正态分布
figure
histfit(x)     %%结果参见图18，加入了较接近的正态分布的密度曲线

图17  例5的hist(x)             图18  例5的histfit(x)
3. 参数估计
当我们可以基本确定数据集符合某种分布时，下一步我们就该确定这个分布的参数了．由于正态分布情况发生的比较多，故一般我们首先考虑的分布将是正态分布．考虑最多的也是正态分布情况．
对于未知参数的估计，可分两种情况：点估计与区间估计．
（1）点估计：构造样本与某个统计量有关的一个函数，作为该统计量的一个估计，称为点估计．Matlab统计工具箱中，一般采用最大似然估计法给出参数的点估计．可以证明：
① 正态分布中，最大似然估计是，的最大似然估计是；
② 泊松分布的最大似然估计是；
③ 指数分布的最大似然估计是，等等．
例6：已知上述例1的数据服从正态分布，试求出和的值．
解：      x=load('data1.txt');
x=x(:);
[mu, sigma] = normfit(x)
mu =
75.3417
sigma =
8.8768
因此，=mu=75.3412，=sigma2=8.87682=78.7982．
（2）区间估计：构造样本与某个统计量有关的两个函数，作为该统计量的下限估计与上限估计，下限与上限一般能够构成一个区间．这个区间作为该统计量的估计，称为区间估计．Matlab统计工具箱中，一般也采用最大似然估计法给出参数的区间估计．
例7：已知上述例1的数据集服从正态分布，试求出和的置信度为95％的区间估计．
解：      x=load('data1.txt');
x=x(:);
[mu, sigma muci, sigmaci] = normfit(x)
mu =
75.3417
sigma =
8.8768
muci =
73.7371
76.9462
sigmaci =
7.8781
10.1678
因此，73.737176.9462，7.878110.1678．
例8：从自动机床加工的同类零件中抽取16件，测得长度值为(data6.txt)：
12.15  12.12  12.01  12.08  12.09  12.16  12.06  12.13
12.07  12.11  12.08  12.01  12.03  12.01  12.03  12.06
已知零件长度服从正态分布，求零件长度的均值和标准差的置信度为99%的置信区间．
解：      x=load('data6.txt');
x=x(:);
[mu, sigma, muci, sigmaci] = normfit(x, 0.01)
mu =
12.0750
sigma =
0.0494
muci =
12.0386
12.1114
sigmaci =
0.0334
0.0892
其中muci(1)、muci(2)分别是平均值在99％置信度下的上下限；而sigmaci(1)、sigmaci(2)分别是标准差在99％置信度下的上下限．
4．正态假设检验
对总体的分布律或分布参数作某种假设，根据抽取的样本观察值，运用数理统计的分析方法，检验这种假设是否正确，从而决定接受假设或拒绝假设，这就是假设检验问题．这里仅以正态假设检验为例，来说明假设检验的基本过程．
正态假设检验的一般过程是：
（1）对比正态分布的概率密度函数图，判断某统计量的分布可能服从正态分布；
（2）利用统计绘图函数normplot( )或weibplot( )进行正态分布检验．
（3）假设检验：利用Matlab统计工具箱给出的常用的假设检验方法的函数ttest(x,m,alpha)，进行显著性水平为alpha的t假设检验，以检验正态分布样本x（标准差未知）的均值是否为m．运行结果中，当h=1时，表示拒绝零假设；当h=0时，表示不能拒绝零假设．
例9：试说明例5所示的刀具的使用寿命服从正态分布，并且说明在方差未知的情况下其均值m取为597是否合理？
解：（1）对比正态分布的概率密度函数图（图17、图18）以及对正态分布的描述（一个变量如果是由大量微小的、独立的随机因素的叠加效果，那么这个变量一定是正态变量．比如测量误差、产品质量等都可用正态分布描述），可得初步结论：该批刀具的使用寿命可能服从正态分布．
（2）利用统计绘图函数normplot(x) 进行分布的正态性检验．由于：
x=load('data5.txt');
x=x(:);
normplot(x)

图19  刀具寿命分布正态性检验
结果如图19所示，经观察这100个离散点非常靠近倾斜直线段，图形为线性的，因此可得出结论：该批刀具的使用寿命近似服从正态分布．
（3）利用函数ttest(x,m,alpha)进行显著性水平为alpha的t假设检验．由于：
x=load('data5.txt');
x=x(:);
h=ttest(x,597,0.05)
得：h = 0
检验结果：h=0，表示不拒绝零假设，说明所提出的假设“寿命均值为597”是合理的．
读者可以验证：当执行h=ttest(x,555,0.05)，将得到h = 1，表示拒绝零假设．请读者自行解释此结果的含义．
四、自己动手
1．了解本实验中虽已提及但没有详细介绍的其余10种概率分布的密度函数，如Beta分布、Gamma分布、Weibull分布等，写出它们的概率分布的密度函数表达式（本实验的附录中已经列出一部分），并画出相应的图形．
2．写出本实验所列出的10种概率累积分布函数表达式，并画出相应的概率累积分布函数图形．
3．用tabulate( )函数将例1、例2的分数数据按频数表的方式进行统计，每5分为一个分数段（可参见例4），观察数据分布有什么规律．
4．用weibplot(x)函数进行例9的正态分布检验，比较与例9的差别．
5．例3给出的上海1998年来的月降雨量的数据(data3.txt) 看上去很接近图10所示的分布，但分布好象没有直接进行参数估计的函数，试寻求对此数据进行参数估计的可能方法．
6．向例3给出的上海1998年来的月降雨量的数据(data3.txt) 中“补充”一些数据，使其看上去很接近正态分布，并求此时的均值和标准差的点估计与置信度为97%的区间估计．
7．在第6题基础上，说明在方差未知的情况下，其均值取为1150是否合理？
8．ttest( )函数的完整用法是：[h,sig,ci] = ttest(x,m,alpha,tail)
其中 sig为观察值的概率，当sig为小概率时则对零假设提出质疑（这里的零假设为：．也可以是其它形式，例如：、等）；ci为真正均值μ的1-alpha置信区间；不写tail，表示其取值为0．
说明：若h=0，表示在显著性水平alpha下，不能拒绝零假设；若h=1，表示在显著性水平alpha下，可以拒绝零假设．
若 tail=0，表示备择（对立）假设为：（默认，双边检验）；若tail=1，表示备择（对立）假设为：（单边检验）；若tail=-1，表示备择（对立）假设为：（单边检验）．
试用该函数求解如下问题：某种电子元件的寿命X（以小时计）服从正态分布，、均未知．现测得16只元件的寿命如下：
159  280  101  212  224  379  179  264  222  362  168  250
149  260  485  170
问当取alpha=0.05时：（1）是否有理由认为元件的平均寿命不大于225（小时）？（2）是否有理由认为元件的平均寿命不大于295（小时）？
9．查看函数 ttest2( )的用法，并用于处理Matlab 统计工具中的数据文件gas.mat．回答问题：一月份油价price1与二月份油价price2的均值是否相同？
五、附录
附录：Matlab中的其它部分概率分布函数名及其数学意义列表：
函数名
对应分布
数学意义
batapdf
Beta分布

chi2pdf
卡方分布

fpdf
F分布


gampdf
Gamma分布
，
raylpdf
瑞利分布
，
tpdf
t分布


weibpdf
Weibull分布