编程手之数学建模学习——数理统计

一、参数估计和检验

1、在一定置信区间上的分布函数参数估计

x0=[506  508  499  503  504  510  497  512
514  505  493  496  506  502  509  496]; x0=x0(:);
pd = fitdist(x0,'Normal')  %对数据进行正态分布拟合
ci = paramci(pd,'Alpha',0.05)  %ci的第一列是均值的置信区间

关于fidist函数中分布名称参数;

 关于常见分布函数的介绍：http://t.csdnimg.cn/znRw2

2、求单侧置信区间上、下限

[h,p,ci]=ttest(x0,xb,'Alpha',0.05,'Tail','right')  %ci是单侧置信区间

3、均值检验、方差检验、总体均值差检验

[h1,p1,ci1,st1]=ttest(x1,mean(x1),'Alpha',0.1) %均值检验和区间估计

[h3,p3,ci3,st3]=vartest(x1,var(x1),'Alpha',0.1) %方差检验和区间估计

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)%总体均值差检验

二、分布假设检验

经验分布函数

当n趋于无穷大时，经验分布函数趋向于总体分布函数。

画出经验分布函数，可以初步判断总体分布是什么；

clc, clear, close all
a=readmatrix('data7_5_1.txt'); a=a(~isnan(a)); %去掉NaN值
[f,x]=ecdf(a) %计算经验分布函数的取值
ecdf(a)  %画经验分布函数图形
grid on, xlabel('$x$','Interpreter','Latex')
ylabel('$F(x)$','Interpreter','Latex')

Q-Q图

Q-Q图可以直观判断假设的拟合优度，多用于正态分布。

clc, clear, close all
a=readmatrix('data7_5_1.txt'); a=a(~isnan(a)); %去掉NaN值
pd=fitdist(a, 'Normal'), qqplot(a,pd)  %Matlab工具箱直接画Q-Q图
sa=sort(a); %把a按照从小到大排列
n=length(a); pi=([1:n]-1/2)/n;
yi=norminv(pi,pd.mu,pd.sigma)'; %计算对应的yi值
hold on, plot(yi,sa,'o') %再重新描点画Q-Q图

卡方拟合检验

用于检验单个总体分布假设的拟合优度，适合非正态分布。

clc, clear, n=100;
bins=[0:7]'; mi=[36  40  19  2   0   2   1   0]'; 
pd = fitdist(bins,'Poisson','Frequency',mi)
expCounts = n * pdf(pd,bins)  %计算期望的频数 %泊松分布
[h,p,st] = chi2gof(bins,'Ctrs',bins,'Frequency',mi, ...
                   'Expected',expCounts,'NParams',1)
k2 = chi2inv(0.95, st.df)  %求临界值,st.df为自由度

pd = 

  PoissonDistribution

  泊松 分布
    lambda = 1   [0.804004, 1.196]

expCounts =

   36.7879
   36.7879
   18.3940
    6.1313
    1.5328
    0.3066
    0.0511
    0.0073

h =

     0

p =

    0.4819

st = 

  包含以下字段的 struct:

    chi2stat: 1.4601
          df: 2
       edges: [-0.5000 0.5000 1.5000 2.5000 7.5000]
           O: [36 40 19 5]
           E: [36.7879 36.7879 18.3940 8.0291]

k2 =

    5.9915

1. pd(泊松分布的参数):
   - `lambda = 1`: 这表明拟合的泊松分布的参数 λ（平均发生率）为 1。
   - 置信区间 [0.804004, 1.196]` 表示我们有 95% 的信心，真实的 λ 值落在这个区间内。

2. **`expCounts` (期望频数)**:
   - 这些是根据拟合得到的泊松分布计算出的期望频数。与观察到的频数 `mi` 进行比较时，这些值用于判断是否存在显著差异。

3. h:
   - `h = 0` 表示不能拒绝零假设。即观察到的频数和期望频数之间没有显著差异。

4. `p`:
   - `p = 0.4819` 表示 p 值。通常，如果 p 值大于 0.05（常用显著性水平），则我们不能拒绝零假设。在这里，p 值为 `0.4819`，远大于 0.05，因此我们认为观察到的频数与期望频数没有显著差异。

5. `st` (统计信息结构体)**:
   - `chi2stat: 1.4601` 是计算出的卡方统计量，这个值用于和临界值进行比较。
   - `df: 2` 是自由度，通常计算方式为 `k - p`，其中 `k` 是分类数，`p` 是估计的参数数量。在此例中，假设有 4 个分类（bin），估计了 1 个参数（λ）。
   - `edges` 是每个 bin 的边界。
   - `O` 是观察到的频数，`E` 是期望频数。

6. `k2`:
   - `k2 = 5.9915` 是在 95% 显著性水平下的卡方临界值。用于和 `chi2stat` 进行比较。

结果解读

根据以上结果，我们可以做出以下结论：

- 由于 `h = 0` 和 `p = 0.4819`，这表明观察到的频数与期望频数之间没有显著差异。因此，我们认为观察数据可以被泊松分布很好地拟合。
- 卡方统计量（`1.4601`）小于临界值（`5.9915`），这进一步支持了这一结论。
- 这表明数据的分布特征与泊松分布是一致的，换句话说，泊松模型是适合此数据的。

 

clc, clear
edges=[0:100:300 inf]; bins=[50 150 250 350]; %定义原始数据区域的边界和中心
mi=[121 78 43 58]; %已知观测频数
pd=makedist('exp',200)  %定义指数分布
expect=sum(mi)*diff(cdf(pd, edges))  %计算期望频数
[h,p,st]=chi2gof(bins,'Edges',edges,'cdf',pd,'Frequency',mi) 
k2=chi2inv(0.95,st.df)  %求临界值，st.df为自由度

秩和检验

用于检验两个总体分布是否相同。取长度短的样本秩和。

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工具箱直接进行检验

三、方差分析

用于比较三个或更多组样本均值是否存在显著差异，常表述为某因素对产品有无显著影响。

单因素方差分析

clc, clear, close all
a=readmatrix('data7_19.txt')' %注意矩阵转置，列
[p,t,st]=anova1(a)
Fa=finv(0.95, t{2,3}, t{3,3})  %计算F分布的上alpha分位数

双因素方差分析

clc, clear, close all
a=readmatrix('data7_20.txt')
[p,t,st]=anova2(a)

四、多元线性回归、多元二项式回归、非线性回归、逐步回归

用于对因素进行拟合并预测。

 多元线性回归：

逐步回归：

总结：

五、灰色预测模型

function []=greymodel(y)
% 本程序主要用来计算根据灰色理论建立的模型的预测值。
% 应用的数学模型是 GM(1,1)。
% 原始数据的处理方法是一次累加法。
y=input('请输入数据 ');
n=length(y);
yy=ones(n,1);
yy(1)=y(1);
for i=2:n
    yy(i)=yy(i-1)+y(i);
end
B=ones(n-1,2);
for i=1:(n-1)
    B(i,1)=-(yy(i)+yy(i+1))/2;
    B(i,2)=1;
end
BT=B';
for j=1:n-1
    YN(j)=y(j+1);
end
YN=YN';
A=inv(BT*B)*BT*YN;
a=A(1);
u=A(2);
t=u/a;
i=1:n+2;
yys(i+1)=(y(1)-t).*exp(-a.*i)+t;
yys(1)=y(1);
for j=n+2:-1:2
    ys(j)=yys(j)-yys(j-1);
end
x=1:n;
xs=2:n+2;
yn=ys(2:n+2);
plot(x,y,'^r',xs,yn,'*-b');
det=0;

sum1=0;
sumpe=0;
for i=1:n
    sumpe=sumpe+y(i);
end
pe=sumpe/n;
for i=1:n;
    sum1=sum1+(y(i)-pe).^2;
end
s1=sqrt(sum1/n);
sumce=0;
for i=2:n
    sumce=sumce+(y(i)-yn(i));
end
ce=sumce/(n-1);
sum2=0;
for i=2:n;
    sum2=sum2+(y(i)-yn(i)-ce).^2;
end
s2=sqrt(sum2/(n-1));
c=(s2)/(s1);
disp(['后验差比值为：',num2str(c)]);
if c<0.35
    disp('系统预测精度好')
else if c<0.5
        disp('系统预测精度合格')
    else if c<0.65
            disp('系统预测精度勉强')
        else
            disp('系统预测精度不合格')
        end
    end
end
            
disp(['下个拟合值为 ',num2str(ys(n+1))]);
disp(['再下个拟合值为',num2str(ys(n+2))]);

[724.57, 746.62, 778.27, 800.8, 827.75,871.1, 912.37, 954.28, 995.01, 1037.2]
%已知量
[2.874,3.278,3.337,3.390,3.679]
%待预测量