MATLAB数理统计------2019/7/24,7/25

数理统计

1. 参数估计和假设检验

• 区间估计
  1） ttest命令是检验某单个样本均值是否等于某个值，进行单个总体，方差未知的 $t$检验，同时给出了参数的区间估计
    用法：
     （1）h = ttest(x) ：零假设为数据的均值$=0$,若数据均值不为０即拒接零假设，则h=1;
      （2）h=ttest(x,m) ：零假设为数据的均值$=m$,若数据均值不为m即拒接零假设，则h=1;
     （3）[h,p,ci,stats] = ttest(x,m,Name,Value) :零假设为数据的均值$=m$,若h=1，则表明在$\alpha$的概率下拒绝零假设;
     Name是参数名称,Value是相应的值。Name必须出现在引号中,可以有多个Name，Value；
    $ci$是置信区间估计
    当Name='Tail’时 若Value='right’则表明检验总体平均值大于m的另一种假设；
  value='left’则表明检验总体平均值小于m的另一种假设；value='booth’则表明检验总体平均值不是m的替代假设；
    stats中：tstat = 测试统计量的值。
         df = 测试的自由度。
        sd = 总体标准差

例7.1
MATLAB命令：tinv(1-alpha/2,n-1) 是 求出$t_{1-\frac{\alpha }{2}}(n-1)$的值
  近似服从正态分布的总体均值$\mu$的一个置信水平为$1-\alpha$的置信区间为$(\overline{X}\underline{+}\frac{S}{\sqrt{n}}{t}_{\frac{\alpha }{2}}(n-1))$
此题中置信度为0.95，所以置信水平为$1-\alpha =0.95$,即$\alpha =0.05$

x0=[506 508 499 503 504 510 497 512 514 505 493 496 506 502 509 496];
x0=x0(:);
alpha=0.05;
mu=mean(x0),sig=std(x0),n=length(x0);
%mean 数组的均值；
%std  算出x的标准偏差
t=[mu-sig/sqrt(n)*tinv(1-alpha/2,n-1),mu+sig/sqrt(n)*tinv(1-alpha/2,n-1)];

[h,p,ci]=ttest(x0,mu,0.05);
% ttest命令是进行单个总体，方差未知的t检验，同时给出了参数的区间估计

例7.2

x0=[1050 1100 1120 1250 1280];
alpha=1-0.95;
mu=mean(x0);  
sig=std(x0);
n=length(x0);
a=mu-sig/sqrt(n)*tinv(1-alpha,n-1)
[h,p,ci]=ttest(x0,mu,alpha)

例7.3

x1=[6.683 6.681 6.676 6.678 6.679 6.672];
x2=[6.661 6.661 6.667 6.667 6.664];
[h1,p1,ci1,st1]=ttest(x1,mean(x1),'Alpha',0.1)
[h2,p2,ci2,st2]=ttest(x2,mean(x2),'Alpha',0.1) 
%均值的检验和区间估计

[h3,p3,ci3,st3]=vartest(x1,var(x1),'Alpha',0.1)  
%方差的检验和区间估计，用法与ttest命令相似ci是置信区间

[h4,p4,ci4,st4]=vartest(x2,var(x2),'Alpha',0.1)

例7.4

x1=[6.683 6.681 6.676 6.678 6.679 6.672];
x2=[6.661 6.661 6.667 6.667 6.664];
[h,p,ci,st]=ttest2(x1,x2,'Alpha',0.1)   
%ttest2用于有两个样本

• 经验分布函数
   1）定义：设$X_1,X_2,\cdots ,X_n$是总体$F$的一个样本，用$S(x)(-\infty \le x\le \infty )$表示$X_1,X_2,\cdots ,X_n$中不大于$x$的随机变量的个数。定义经验分布函数$F_n(x)$为
  $$F_n(x)=\frac{1}{n}S(x),-\infty \le x\le \infty$$
   2）$F_n(x)$的观察值：（$F_n(x)$的观察值仍以$F_n(x)$表示）
    一般地，设$x_1,x_2,\cdots ,x_n$是总体$F$的一个容量为$n$的样本值。先将$x_1,x_2,\cdots ,x_n$按自小到大的次序重新排列编号。设为：
  $$x_{(1)}\le x_{(2)}\le \cdots \le x_{(n)}$$
  则经验分布函数$F_n(x)$的观察值为：
  $$F_n(x)=\{\begin{array}{l}0,x<x_{(1)}\\ \frac{k}{n},若x_{(k)}\le x\le x_{(k+1)}\\ 1,若x\ge x_{(n)}\end{array}$$
  例7.5（经验分布函数的应用）

a=[141 148 132 138 154 142 150 146 155 158
   150 140 147 148 144 150 149 145 149 158
   143 141 144 144 126 140 144 142 141 140
   145 135 147 146 141 136 140 146 142 137
   148 154 137 139 143 140 131 143 141 149
   148 135 148 152 143 144 141 143 147 146
   150 132 142 142 143 153 149 146 149 138
   142 149 142 137 134 144 146 147 140 142
   140 137 152 145  0   0    0  0   0   0]
a=nonzeros(a);
%去除多余的0，并将矩阵a展成列向量

[ycdf,xcdf,n]=cdfcalc(a) 
%计算经验分布函数的取值 ,ycdf为当x==xcdf时的概率，xcdf为x从小到大的排序（无重复），
%n为（含重复的x的取值）一共的x数,比如此题中a中有84个值，所以n=84

cdfplot(a),title('')  

%画CDF曲线，CDF指在x点左侧事件发生的总和，即在曲线上一点（x1,y1）处标明x<=x1的概率为y1

hold on,plot(xcdf,ycdf(2:end),'.')
%再重新画经验分布函数的取值

• Q-Q图
   Q-Q图用来检验拟合度
   1）基本思想：将经验分布函数的分位数点和分布模型的理论分位数点作为一对数组画在坐标图上，如果这$n$个点看起来像一条直线，说明拟合效果好。
   2）Q-Q图计算步骤：
    (1)将$x_1,x_2,\cdots ,x_n$按自小到大的次序重新排列为：
  $x_{(1)}\le x_{(2)}\le \cdots \le x_{(n)}$
    (2)取$y_i=F^{-1}((i-1/2)/n),i=1,2,\cdots ,n$
    (3)将$(y_i,x_{(i)}),i=1,2,\cdots ,n$这$n$个点画为直角坐标系上
    (4)若这$n$个点看起来成一条$45°$角的直线，则$x_1,x_2,\cdots ,x_n$拟合分布$F(x)$效果好
  例7.6

a=[141 148 132 138 154 142 150 146 155 158
   150 140 147 148 144 150 149 145 149 158
   143 141 144 144 126 140 144 142 141 140
   145 135 147 146 141 136 140 146 142 137
   148 154 137 139 143 140 131 143 141 149
   148 135 148 152 143 144 141 143 147 146
   150 132 142 142 143 153 149 146 149 138
   142 149 142 137 134 144 146 147 140 142
   140 137 152 145  0   0    0  0   0   0]
a=nonzeros(a);
xbar=mean(a);
%求均值

s=std(a);
%求标准差

pd=makedist('normal','mu',xbar,'sigma',s)
%利用makedist命令定义正态分布
%makedist('distname','name','value')表示创建一个概率分布函数：
%当'distname'=='normal'时表示创建一个正太分布函数，参数名为'name',参数值为'value',
%可以有多个参数名和参数值


qqplot(a,pd)  %使用MATLAB工具箱直接画Q-Q图
%qqplot(x,pd)显示样本数据x的分位数的分位数图与指定分布的理论分位数概率分布对象pd。
%如果x的分布与pd指定的分布相同,则绘图显示为线性。



%不利用工具箱画Q-Q图：
sa=sort(a);
n=length(a);
pi=([1:n]-1/2)/n;
yi=norminv(pi,xbar,s)';
%x=norminv(p,mu,singa)命令求正态分布分位数：若p{X≥xa}=α,则称X的上侧α分位数为xa

hold on,plot(yi,sa,'.')

• 非参数检验
  1）${\chi }^2$拟合优度检验：
   步骤：
   (1)建立待检验假设$H_0:$总体$X$的概率密度为$f(x)$
   (2)在数轴上选取$k-1$个分点$t_1,t_2,\cdots ,t_{k-1}$,将数轴分为$k$个区间$:(-\infty ,t_1),[t_1,t_2),\cdots ,[t_{k-2},t_{k-1}),[t_{k-1},+\infty ),$令$p_i$为分布函数$F(x)$的总体$X$在第$i$个区间内取值的概率，设$m_i$为$n$个样本观察值中落入第$i$个区间上的个数，也称组频数。
   (3)选取统计量${\chi }^2={\sum }_{i=1}^k\frac{(m_i-n{p}_i{)}^2}{n{p}_i}$,若$H_0$为真，则${\chi }^2\sim {\chi }^2(k-1-r)$,r为分布函数$F(x)$中未知参数的个数
   (4)对于给定的显著水平$\alpha$,确定${\chi }_{\alpha }^2$使其满足 P { χ 2 ( k − 1 − r ) &gt; χ 2 } = α P\{\chi^{2}(k-1-r)&gt;\chi^{2}\}=\alpha P{χ2(k−1−r)>χ2}=α,并依据样本计算统计量${\chi }^2$的观察值。
   (5)做出判断：若${\chi }^2<{\chi }_{\alpha }^2,$则接受$H_0$
  例7.7（判断样本是否服从泊松分布）

n=100;
f=0:7;
num=[36 40 19 2 0 2 1 0];
lamda=dot(f,num)/100;;%dot 为计算点积,
%此题中f,num为实数向量故计算结果为f(1)*num(1)+f(2)*num(2)+...+f(end)*num(end)
%若f,num为矩阵则将f和num的各列视为向量，并计算对应列的点积


pi=poisspdf(f,lamda);%计算参数为lamda时泊松分布在f处的值
[h,p,st]=chi2gof(f,'ctrs',f,'frequency',num,'expected',n*pi,'nparams',1)
%[h,p,st]=chi2gof(x,name1,val1,name2,val2,.......)
col3=st.E/sum(st.O); 
%计算分布函数F(x)的总体X在第i个区间内取值的概率


col4=st.E;  
%显示:每个小区间对应的理论频数

col5=st.O.^2./col4;   
%显示:实际落入一个区间的个数的平方/该区间的理论频数。
%st.O:落入每个小区间内观测的个数，即时间频数                                                                  

sumcol5=sum(col5)  
%col5求和，χ2=sumcol5的个位数及其之后的数.eg:若sumcol5=101.4611,则χ2=11.4611
k2=chi2inv(0.95,st.df)

%求χ2(k−1−r)的值,若χ2(k−1−r)>χ2,则接受假设H0，即样本来自泊松分布的总体

MATLAB函数：

 [ h,p,st ]=chi2gof ( x,name1,val1,name2,val2,…)
  h: 等于0或1，若h=0则在显著性水平α下接受原假设
 p: 当p值小于或等于显著性水平α时，拒绝原假设，否则接受原假设。
 st:包含以下字段：

chi2stat:卡方检验统计量
df ： 自由度
edges：合并后各区间的边界向量
O： 落入每个小区间内观测的个数，即时间频数
E： 每个小区间对应的理论频数

参数名	参数值	说明
‘nbins’	正整数，默认值为10	分组（或区间）的个数
‘edges’	向量	指定各个区间的边界
‘ctrs’	向量	指定各区间的中点
上述三个参数不能同时指定，一次调用最多只能指定其中的一个参数，因为后面两个参数已经潜在的指定了分组数了
cdf	概率分布对象 / 功能手柄 / 单元格数组	如果CDF是概率分布对象,则自由度将考虑是使用fitdist估计参数还是使用makedist参数. 如果CDF是函数句柄,则分步函数必须将x作为其唯一的参数。 如果CDF是单元格数组,则第一个元素必须是函数句柄,其余元素必须是参数值,每个单元格一个。函数必须将x作为其第一个输入参数,将数组中的其他参数作为后序的输入参数。
expected’	向量	指定各区间的理论频数，与‘cdf’不能同时出现
NParams’	向量	指定分布中待估参数的个数，它确定了卡方分布的自由度
emin’	非负整数，默认值5	指定一个区间对应的最小理论频数，初始分组中， 理论频数小于这个值的区间和相邻区间合并。如果指定为0， 将不进行区间合并
‘frequency’	与x等长的向量	指定x中各元素出现的频数
‘alpha’	0--1之间的数，默认值0.05	指定

例7.8

edges=[0:100:300 inf];
bins=[50 150 250 inf];
num=[127 78 43 58];
pd=makedist('exp',200); 
%定义指数分布

[h,p,st]=chi2gof(bins,'Edges',edges,'cdf',pd,'Frequency',num)
pi=st.E/sum(st.O);
co14=st.E;  
%显示:每个小区间对应的理论频数

col5=st.O.^2./col4   
%显示:实际落入一个区间的个数的平方/该区间的理论频数。st.O:落入每个小区间内观测的个数，即时间频数

sumcol5=sum(col5)
k2=chi2inv(0.95,st.df)
%求χ2(k−1−r)的值,若χ2(k−1−r)>χ2,则接受假设H0，即样本来自泊松分布的总体