第七章：数理统计

最新推荐文章于 2021-10-10 19:44:47 发布

victor_cs_

最新推荐文章于 2021-10-10 19:44:47 发布

阅读量1.5k

点赞数 3

分类专栏：数学建模算法与应用文章标签：数理统计数学建模

本文链接：https://blog.csdn.net/qq_41000485/article/details/99423456

版权

数学建模算法与应用专栏收录该内容

16 篇文章 25 订阅

订阅专栏

7.1 参数估计和假设检验

7.1.1 区间估计

clc, clear
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);
t=[mu-sig/sqrt(n)*tinv(1-alpha/2,n-1),mu+sig/sqrt(n)*tinv(1-alpha/2,n-1)]
%以下命令ttest的返回值ci就直接给出了置信区间估计
[h,p,ci]=ttest(x0,mu,0.05)  %通过假设检验也可求得置信区间

ttest命令即进行单个总体，方差未知的t检验，同时给出了参数的区间估计。

clc, clear
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) %方差检验和区间估计
[h4,p4,ci4,st4]=vartest(x2,var(x1),'Alpha',0.1)

7.1.2 经验分布函数

对于任一实数x，当n充分大时，经验分布函数的任一个观察值与总体分布函数只有微小的差别，从而在实际上可当做总体分布函数来使用。

clc, clear
a=textread('ex7_5.txt'); a=nonzeros(a); %读入数据，并去掉多余的零展成列向量
[ycdf,xcdf,n]=cdfcalc(a) %计算经验分布函数的取值
cdfplot(a), title('') %画经验分布函数的图形
hold on, plot(xcdf,ycdf(2:end),'.') %再重新画经验分布函数的取值
xlswrite('ex7_5.xls',[xcdf,ycdf(2:end)])

7.1.3 Q-Q图

Q-Q图是检验拟合优度的，如果看起来成一条45度角的直线，则拟合分布的效果很好。

clc, clear
a=textread('ex7_5.txt'); a=nonzeros(a); 
xbar=mean(a), s=std(a) %求均值和标准差
pd=ProbDistUnivParam('normal',[xbar s]) %定义正态分布
qqplot(a,pd)  %Matlab工具箱直接画Q-Q图
%下面不利用工具箱
sa=sort(a); %把a从小到大排列
n=length(a); pi=([1:n]-1/2)/n;
yi=norminv(pi,xbar,s)' %计算对应的yi值
hold on, plot(yi,sa,'.') %重新描点画图

7.1.4 非参数检验

（1）卡方拟合优度检验

clc, clear, n=100;
f=0:7; num=[36  40  19  2   0   2   1   0]; 
lamda=dot(f,num)/100
pi=poisspdf(f,lamda) 
[h,p,st]=chi2gof(f,'ctrs',f,'frequency',num,'expected',n*pi,'nparams',1) %调用工具箱
col3=st.E/sum(st.O) %计算表中的第3列数据
col4=st.E %计算表中的第4列数据
col5=st.O.^2./col4  %计算表中的第5列数据
sumcol5=sum(col5)  %计算表中的第5列数据的和
k2=chi2inv(0.95,st.df)  %求临界值，st.df为自由度

（2）柯尔莫哥洛夫检验

clc, clear
a=textread('ex7_5.txt'); a=nonzeros(a); 
xbar=mean(a), s=std(a), s2=var(a) %
[yn,xn]=cdfcalc(a); 
yn(end)=[];  %yn的元素个数比xn多了一个，删除最后一个值
y=normcdf(xn,xbar,s); %计算理论分布函数值
Dn=max(abs(yn-y))  %计算统计量的值
LJ=1.36/sqrt(length(a)) %计算拒绝域的临界值
%以下直接调用Matlab工具箱的命令进行KS检验
pd=makedist('Normal','mu',xbar,'sigma',s)
[h,p,st]=kstest(a,'CDF',pd) %直接调用工具箱的命令进行KS检验

7.1.5 秩和检验

检验两个总体有相同的分布

clc, clear
x=[41.5	42.0	40.0	42.5	42.0	42.2	42.7	42.1	41.4];
y=[41.2	41.8	42.4	41.6	41.7	41.3];
yx=[y,x]; yxr=tiedrank(yx) %计算秩
yr=sum(yxr(1:length(y))) %计算y的秩和
[p,h,s]=ranksum(y,x) %利用Matlab工具箱直接进行检验

7.2 Bootstrap方法

7.2.1 非参数

设总体的分布F未知，但已知有一个容量为n的来自分布F的数据样本，这一样本按放回抽样的方法抽取一个容量为n的样本，这个样本即Bootstrap样本或资助样本。
（1）估计量的标准误差的Bootstrap估计

clc, clear
a=[9.5  18.2  12.0  10.2  18.2
21.1  18.2  12.0  9.5  10.2
21.1  10.2  10.2  12.0  10.2
18.2  12.0  9.5  18.2  10.2
21.1  12.0  18.2  12.0  18.2
10.2  10.2  9.5  21.1  10.2
9.5  21.1  12.0  10.2  12.0
10.2  18.2  10.2  21.1  21.1
10.2  10.2  18.2  18.2  18.2
18.2  10.2  18.2  10.2  10.2];
zw=quantile(a',0.5) %求矩阵a每一行的中位数
bzc=std(zw)  %计算中位数标准差的Bootstrap估计

clc, clear
a=[18.2  9.5  12.0  21.1  10.2]; %输入原始样本
b=bootstrp(1000,@(x)quantile(x,0.5),a) %计算各个Bootstrap样本的中位数
c=std(b)  %计算中位数标准差

（2）估计量的均方误差的Bootstrap估计

clc, clear
a=[136.3  136.6  135.8  135.4  134.7  135.0  134.1  143.3  147.8  148.8  134.8  135.2  134.9  149.5  141.2  135.4  134.8  135.8  135.0  133.7  134.4  134.9  134.8  134.5  134.3  135.2];
b=bootstrp(10000,@(x)quantile(x,0.5),a);  %计算各个Bootstrap样本的中位数
c=mean((b-quantile(a,0.5)).^2)  %求均方误差

（3）Bootstrap置信区间

clc, clear
a=[9  8  10  12  11  12  7  9  11  8  9  7  7  8  9  7  9  9  10  9  9  9  12  10  10  9  13  11  13  9];
b=bootci(10000,{@(x)[mean(x),std(x)],a},'alpha',0.1) 
%返回值b第一列为均值的置信区间，第二列为标准差的置信区间

7.2.2 参数

（1）求置信下限

clc, clear
a=[142.84  97.04  32.46  69.14  85.67  114.43  41.76  163.07  108.22  63.28];
eta=sqrt(mean(a.^2))  %求最大似然估计
beta=2; B=5000; alpha=0.05;
b=wblrnd(eta,beta,[B,10]);  %产生服从韦布尔分布的随机数
etahat=sqrt(mean(b.^2,2));  %计算每个样本对应的最大似然估计
seta=sort(etahat);  %把B个最大似然估计按照从小到大排列
k=floor(B*alpha)
se=seta(k)  %提取相应位置的估计量
Rt0=exp(-(50/se)^2)  %求对应的置信下限

（2）求置信区间

clc, clear
x0=[342  500  187];
theta=(x0(2)+2*x0(3))/sum(x0)/2  %计算最大似然估计
fb=[(1-theta)^2,2*theta*(1-theta),theta^2]  %计算分布律
cf=cumsum(fb)  %求累计分布
a=rand(1029,1000);  %每一列随机数对应一个bootstrap样本
jx1=(a<=cf(1));  %1对应M出现
jx2=(a>cf(1) & a<=cf(2)); %1对应MN出现
jx3=(a>=cf(2)); %1对应N出现
x1=sum(jx1); x2=sum(jx2); x3=sum(jx3);  
theta2=(x2+2*x3)/1029/2;  %计算统计量theta的值
stheta=sort(theta2); %把统计量按照从小到大排序
qj=[stheta(50), stheta(950)]  %提出置信区间的取值

7.3 方差分析

单因素试验的方差分析如下，双（多）因素试验的方差分析是类似的

clc, clear, alpha=0.05;
a=[41	65	45
48	57	51
41	54	56
49	72	48
57	64	48];
[p,t,st]=anova1(a) %返回值t是细胞数组
F=t{2,5}  %显示F统计量的值
fa=finv(1-alpha,t{2,3},t{3,3})  %计算临界值

7.4 回归分析

7.4.1 多元线性回归

7.4.2 多元二项式回归

统计工具箱提供了一个rstool的命令
rstool(X,Y,model,alpha)
alpha为显著性水平
model为

linear 线性
purequadratic 纯二次
interaction 交叉
quadratic 完全二次

clc, clear
ab=textread('ex7_19.txt'); 
y=ab(:,[2:5:10]); %提取因变量y的观察值
Y=nonzeros(y); %去掉y后面的0，并变成列向量
x123=[ab([1:13],[3:5]); ab([1:12],[8:10])]; %提取x1,x2,x3的观察值
X=[ones(25,1),x123]; %构造多元线性回归分析的数据矩阵X
[beta,betaint,r,rint,st]=regress(Y,X)  %计算回归系数和统计量等，st的第2个分量就是F统计量，下面根据统计量的表达式重新计算的结果和这里是一样的
q=sum(r.^2)   %计算残差平方和
ybar=mean(Y)  %计算y的观察值的平均值
yhat=X*beta;   %计算y的估计值
u=sum((yhat-ybar).^2)   %计算回归平方和
m=3;    %变量的个数，拟合参数的个数为m+1
n=length(Y); %样本点的个数
F=u/m/(q/(n-m-1))  %计算F统计量的值，自由度为样本点的个数减拟合参数的的个数
fw1=finv(0.025,m,n-m-1) %计算上1-alpha/2分位数
fw2=finv(0.975,m,n-m-1) %计算上alpha/2分位数
c=diag(inv(X'*X))   %计算c(j,j)的值
t=beta./sqrt(c)/sqrt(q/(n-m-1))  %计算t统计量的值
tfw=tinv(0.975,n-m-1)   %计算t分布的上alpha/分位数
save xydata Y x123  %把Y和x123保存到mat文件xydata中，供问题(3)的二次模型使用

clc, clear
load xydata
rstool(x123,Y)

7.4.3 非线性回归

使用统计工具箱中的nlinfit，nlparci，nlpredci，nlintool命令

clc, clear
xy0=[8.55	470	300	10
3.79	285	80	10
4.82	470	300	120
0.02	470	80	120
2.75	470	80	10
14.39	100	190	10
2.54	100	80	65
4.35	470	190	65
13.00	100	300	54
8.50	100	300	120
0.05	100	80	120
11.32	285	300	10
3.13	285	190	120];
x=xy0(:,[2:4]);
y=xy0(:,1);
huaxue=@(beta,x) (beta(4)*x(:,2)-x(:,3)/beta(5))./(1+beta(1)*x(:,1)+...
beta(2)*x(:,2)+beta(3)*x(:,3)); %用匿名函数定义要拟合的函数
beta0=[0.1,0.05,0.02,1,2]';  %回归系数的初值，可以任意取，这里是给定的
[beta,r,j]=nlinfit(x,y,huaxue,beta0)  %计算回归系数beta;t,j是下面命令用的信息
betaci=nlparci(beta,r,'jacobian',j)  %计算回归系数的置信区间
[yhat,delta]=nlpredci(huaxue,x,beta,r,'jacobian',j) %计算y的预测值及置信区间半径

也可以用nlintool得到交互式画面，用Export传送数据
使用命令nlintool(x,y,huaxue,beta0)

7.5 基于灰色模型和Bootstrap理论的大规模定制质量控制方法研究

灰色模型GM(1,1)，适用于大规模定制生产的质量预测分析。

所需信息较少
不需要原始数据分布的先验特征
可保持原系统特征

通过Bootstrap方法，重复抽样，获得一定规模的样本量，进而得到统计量的经验分布并进行区间估计

clc, clear
x0=[2.5320,2.6470,2.6290,2.5840,2.6090,2.6010,2.5280,2.5630,2.6540,2.6190];
n=length(x0);
me=quantile(x0,0.5)  %计算中位数
[h,p,stat]=runstest(x0,me)  %进行游程检验
x1=cumsum(x0);  %求累加序列
zk=(x1(1:end-1)+x1(2:end))/2  %求累加序列的均值序列
B=[-zk', ones(size(zk'))]; yn=x0(2:end)';
ab=B\yn  %拟合参数a，b
syms x(t)
x=dsolve(diff(x)+ab(1)*x==ab(2),x(0)==x0(1)); %求微分方程的符号解
xx=vpa(x,6)  %显示小数格式的符号解
yuce=subs(x,'t',[0:n+5]);  %求累计序列的预测值
yuce=double(yuce); %吧符号数转换成数值类型
yuce0=[x0(1),diff(yuce)]  %求原始数据的预测值
c=std(yuce0(1:n))/std(x0)  %求后验差比值c
nyuce=yuce0(n+1:end)  %提取6个新的预测值
nyb=[x0, nyuce];  %构造新的样本数据
nnyb=reshape(nyb,[4,4])
mu=mean(nnyb)  %分别求4个子样本的均值
jc=range(nnyb)  %分别求4个子样本的极差
xlswrite('hb.xls',[nnyb;mu;jc])  %把数据写到Excel文件中，便于做表使用
b=rand(4,1000); %产生4行1000列的随机数矩阵
h=floor(b*length(nyb))+1; %把随机数映射为编号（每列对应Bootstrap样本编号）
bb=repmat(nyb',1,1000); bb=bb(h); %对新序列进行重复抽样
mmu=mean(bb); mjc=range(bb); %计算1000个子样本的均值和极差
smu=sort(mmu); sjc=sort(mjc); %把均值和极差按照从小到大的次序排列
alpha=0.0027; k1=floor(1000*alpha/2), k2=floor(1000*(1-alpha/2))
mqj=[smu(k1), smu(k2)]  %显示均值的置信区间
jqj=[sjc(k1), sjc(k2)]  %显示极差的置信区间
subplot(1,2,1), plot(mu,'*-'), hold on, plot([1,4],[mqj(1),mqj(1)])
plot([1,4],[mqj(2),mqj(2)]), ylabel('样本均值')
subplot(1,2,2), plot(jc,'*-'), hold on, plot([1,4],[jqj(1),jqj(1)]), 
plot([1,4],[jqj(2),jqj(2)]), ylabel('极差')