【数学建模】基于灰色模型和Bootstrap的集成方法的应用

模型优势

（1）灰色模型在极小样本量情况下进行质量数据预测具有独特的优势，预测效果也相对较好；
（2）基于Bootstrap理论的统计推断方法通过重复抽样，能对未知分布的随机变量的分布参数进行较为精确的区间估计，为构建质量控制图提供依据。

案例分析

某航空产品制造厂生产的一批框段根据技术指标及安装位置的差异性，其半径、形状、连接方式各不相同，共有20种，且每一种产品批生产数量都不超过30件，其生产方式可视为大规模定制模式。

分析流程：

游程检验

游程检验亦称“连贯检验”, 是根据样本容量和游程的多少来判定一个给定样本序列是否排列随机的检验方法。
在以一定顺序（如时间）排列的有序数列中，具有相同属性（如符号）的连续部分被称为一个游程，一个游程中所包含数据的个数称为游程的长度，通常用R表示一个数列中的游程总数。
采用该方法对初始样本数据的随机性进行检验。若排列随机，表明初始样本对应的生产过程状态正常，进行质量数据预测有意义。若排列非随机，应根据5M1E分析法及时检查并纠正系统性质量因素。

5M1E分析造成产品质量的波动的主要原因6个因素的总称：
(1) 人(Man/Manpower) 操作者对质量的认识、技术熟练程度、身体状况等；
(2) 机器(Machine) 机器设备、工夹具的精度和维护保养状况等；
(3) 材料(Material) 材料的成分、物理性能和化学性能等；
(4) 方法(Method) 这里包括加工工艺、工装选择、操作规程等；
(5) 测量(Measurement) 测量时采取的方法是否标准、正确；
(6) 环境(Environment) 工作地的温度、湿度、照明和清洁条件等；

matlab里自带的函数runstest进行游程检测,格式：

[h,p,stats]=runstest(x,v);
[h,p,stats]=runstest(x,′ud′);
[h,p,stats]=runstest(x,′alpha′);

其中：v为临界值,基于它进行游程检验。默认值为数据序列x的均值; 选项′ud′根据数列中出现的连续增和连续减的子序列数确定游程总数,从而进行游程检验。此时数列x中与前一元素相同的数将会被去除;返回值中h为判断结果;p为接受概率,当p值小于或等于显著性水平ahpha时,拒绝原假设,否则接受原假设;stats分4个部分：nruns为游程数；n1为取值1的数目(数据序列中大于v的数据个数)；n0为取值0的数目；Z为检验值。

游程检验临界值表


matlab程序

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);

运行结果

me =2.6050；
h =0；p =1
stat = 包含以下字段的struct:
nruns: 6 %游程数
n1: 5 %取值1总数
n0: 5%取值0总数
z: 0 %检验值

灰色模型求解白化方程参数并求预测值

以新的样本序列中数据为原始数据, 计算样本均值和样本极差的控制限。依时间顺序将其分为 4 个子样本组。
在 Matlab 中用 Bootstrap 工具箱，基于样本数列进行重复抽样，子样本容量为4，抽取1000 个子样本，计算子样本均值和极差。
将得到的1000 个子样本的均值和极差按从小到大的顺序排列，则得到Bootstrap 经验分布。

根据样本均值和样本极差控制图上下限和样本统计量取值，分别绘制样本均值控制图（左）和样本极差控制图（右），观察样本统计量。

数据点均在控制界限内，可以判断生产过程受控。

完整代码

x0=[2.5320,2.6470,2.6290,2.5840,2.6090,2.6010,2.5280,2.5630,2.654
0,2.6190];
n=length(x0);
me=quantile(x0,0.5);
[h,p,stat]=runstest(x0,me);
x1=cumsum(x0); %生成AGO
zk=(x1(1:end-1)+x1(2:end))/2; %生成均值序列z(k)
B=[-zk',ones(size(zk'))]; %生成矩阵B(7.26)
yn=x0(2:end)'; %生成向量序列Y(7.26)
ab=B\yn; %由（7.25）求拟合参数a,b
syms x(t)
x=dsolve(diff(x)+ab(1)*x==ab(2),x(0)==x0(1)); %求解方程(7.24)
yuce=subs(x,'t',[0:n+5]); %求累加数据预测值
yuce=double(yuce);
yuce0=[x0(1),diff(yuce)]; %求原始数据预测值
c=std(yuce0(1:n))/std(x0); %求后验差比值c(7.27)
nyuce=yuce0(n+1:end);
nyb=[x0,nyuce]; %构造新样本, 容量为16
nnyb=reshape(nyb,[4,4]); %将nyb依列重新排成4*4形式
mu=mean(nnyb); %依列求均值
jc=range(nnyb); %依列求极差(最大值与最小值之差)
xlswrite('例7_5_4',[nnyb;mu;jc]);
b=rand(4,1000);
h=floor(b*length(nyb))+1; %通过取整,映射为1到16的编号
bb=repmat(nyb',1,1000); %将nyb转置,重复取1000列
bb=bb(h); %按照h的编号，对bb进行抽样，得到1000个新样本
mmu=mean(bb);mjc=range(bb);
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)])
subplot(1,2,2),plot(jc,'*-'),hold on,plot([1,4],[jqj(1),jqj(1)])