横向最小二乘拟合与纵向卡尔曼滤波方法文档格式.docx
《横向最小二乘拟合与纵向卡尔曼滤波方法文档格式.docx》由会员分享，可在线阅读，更多相关《横向最小二乘拟合与纵向卡尔曼滤波方法文档格式.docx（10页珍藏版）》请在冰豆网上搜索。
02025.995.995.395.395.395.394.194.193.592.992.41.801670
02034.433.83.83.83.83.83.162.532.531.270.63001580
020413.7712.1511.619.729.187.836.755.674.322.431.350.5403704
020536.7834.6831.5329.4327.0625.2223.1221.8118.1316.5513.48.933.943806
020641.5839.1836.0832.9931.6228.8724.7423.0218.915.4613.49.284.472910
0207161472.4969.3962.5854.5247.7143.9940.2734.730.3626.6422.313.013.720208198575.5771.5469.0264.4856.9352.945.3436.7828.2120.6513.67.561.5102092671112.32110.45108.57104.0895.8484.6174.8865.8952.0442.3127.3311.231.8702102107121.97119.6116.28115.33107.7496.3584.4869.2954.1139.3922.7811.392.850211139995.7895.7894.3592.2185.7882.272.1961.4747.1840.0325.7312.873.570212403101.74101.7494.2991.8189.3384.3781.896752.1144.6732.267.447.4403016450122.79122.79122.48121.55119.84115.5108.0698.2982.6466.9844.9622.023.7203022522143.93143.93143.93143.93141.95139.57135.21125.69106.6684.4662.2525.381.590303290060.3460.3460.3460.346058.2855.8651.7246.2133.116.551.03
0304112718.6318.6318.6318.6318.6316.8615.9713.317.992.660
030581814.6714.6714.6714.6713.4513.4513.45118.561.22
030611995.845.845.845.845.845.8451.670
0307183113.6513.6513.6513.6513.1110.387.10.55
030817545.75.75.75.74.561.710
030921630.920.920.920.920.460.46
0310238900000
031124340000
03121171000
提示:
1.预测时用的数据表最好是增量表，就是把原表相邻列作差的到的表，含义是第几个月期间的千车故障数。
预测后再恢复到原表的形式。
2.轿车出厂后的运输是个复杂的事，体积大又贵重，要花费很多时间，从表中数据分析可以得到:
出厂后三个月才开始有销售量，于是每个批次的前三个数据（斜三列）可认为是无效数据。
2
采用横向最小二乘拟合与纵向卡尔曼滤波方法的联合预测方法对原始数据处理进行处理:
数据处理步骤如下:
步骤1:
基于分析结果2,出厂后三个月才有销售量.去除原始表中的斜三列中得数据.结果如下:
表二去除斜三列后数据表（节选）
1211109876543210
„„„„„„„„„„„„
021195.7895.7894.3592.2185.7882.272.1961.4747.1840.0325.7312.873.570212101.7494.2991.8189.3384.3781.896752.1144.6732.267.447.440301122.48121.55119.84115.5108.0698.2982.6466.9844.9622.023.720302143.93141.95139.57135.21125.69106.6684.4662.2525.381.59030360.346058.2855.8651.7246.2133.116.551.03030418.6318.6316.8615.9713.317.992.660030514.6713.4513.4513.45118.561.2203065.845.845.8451.670030713.6513.1110.387.10.5503085.74.561.71003090.920.460.46031000031100312
步骤2:
表的修正
故障部件数，1000
（1）原始的千车故障数=1迄今为止的汽车售出量
故障部件数，1000修正后的千车故障数=2满足使用月数条件的售出量
3
以0205批次使用月数为10个月解释式2分母的“满足使用月数条件的售出量”,0205批次的汽车要到2002年09月份才有销售量,而在2003年6月份以后的售出量（不包括该月）到2004年03月份为止的使用月数还不到10个月,因此满足月数条件的售出量是2002年09月到2003年06月份的销售量。
（2）比较上面1和2式,发现两式的的共同之处在于有一样的“故障部件数”,又基于每月销售量相同的假设下,不难得出由原始的千车故障数向修正的千车故障数的转化与总销售量无关,仅仅与月数有关，关系如下
该批次到制表时的总销售月数修正后的千车故障数=原始的千车故障数，该批次到制表时满足条件的月数（3）修正算法
总销售月数为24-;
（前三个月没有销售量）i
满足条件的月数25--。
ji
修正算法如下:
For=1:
24i
13j
If+<
25ji
（24,i）=，//注意:
使用月数-1不是//，，，，jjszysi,jysi,j（25,i,j）
Else
=0//表格中的空数据赋为0//，，szysi,j
由该算法得到的修正数据见附表1。
步骤3:
对表中的列作残差,也即把原来相邻的列作差得到新的增量表,表示第几个月期间的千车故障数。
步骤4:
基于分析结果1,去掉=13的行.（见附表2）至此在以后的预测计算i
中,0301批次以后的数据都向上挪一行,例如，预测0306批次时,=17,预测0310i批次时=21。
i
步骤5:
表内数据的横向最小二乘拟合与纵向卡尔曼（kalman）滤波方法的联合预测
对于修正后的差分表,同一个批次在相邻月份内发生的千车故障率必然有相关系数,而且故障率可以认为是线性关系,因此横向采取线性最小二乘拟合未
4
知的故障数,再在此基础上运用纵向kalman滤波对拟合后的数据进行除噪处理,从而降低了数据的误差。
例如，对于（0212,13）未知“故障数”用（0212,0）（0212,1）„..（0212,12）数据线性最小二乘法拟合得到，然后通过对（0201,13）,（0202,13）„（0211,13）,（0212,13）进行kalaman滤波分析修正最小二乘法拟合得到的（0212,13）值。
表三横向最小二乘拟合与纵向卡尔曼滤波方法的联合预测顺序表
销售生产制表
月份月份时销
售量
030302111399622.57213.72109.1776.65444.92844.03334.17831.7416.1121.63116.46610.3723.5703040212403655.14198.599.2565.50238.70848.92336.69218.60520.84830.5960.676367.44（11）
030603022522729.55244.52127.2186.64573.61357.1142.84549.61226.611.59（12）
（1）
030703032900273.0695.1649.15532.58923.78126.75823.93817.5891.03（13）
（2）
03080304112774.5229.5613.0210.64410.6437.61333.040（14）（3）
0309030581855.61515.6927.84588.13755.41338.76671.22（15）（4）？
03100306119917.525.844.185.4962.0040（16）（5）？
03110307183135.47515.4758.4258.3250.55（17）（6）？
03120308175413.686.842.280（18）（7）0401030921632.070.230.46（19）（8）04020310238900（20）（9）0403031124340（21）（10）040403121171（22）（11）
具体处理过程如下:
（1）从空表项的最上的一条对角线开始。
用最小二乘法拟合0302批次使用月数为10的数据（=13,=11）ji
（2）用纵向滤波对（=13,=11）的数据进行除噪处理,得到修正值.ji
（3）=+1，=-1，重复进行
（1），
（2）的步骤。
直到该对角线填满为止。
jjii
（4）对下面的对角线，重复进行
（1），
（2），（3）的步骤直到，表中的空表项填满为止。
（其数据处理的顺序如表6.1中的数字所示）
至此,数据处理部分全部结束,得到的数据表中的数据称之为对应批次对应月的“故障数”，以下的数据建模和预测都是基于“故障数”的基础上。
附录一:
5
KALMAN滤波原理
滤波是从获得的测量信息中尽可能的滤除干扰，分离所需要的真实信号。
例
如:
已知测量向量序列Y，Y，„.,，如何求得向量X的估计。
由于燥声的干12k扰，不可能精确算得状态向量X的真值，而只能在一定的统计准则下作出最优k
估计，kalman滤波采用的是最小均方误差估计原则记:
ttttY（k）=（Y,Y,„,Y）12k
~^
希望由Y（k）对j时刻的状态X进行估计。
记估计量和估计误差分别是XjXjj?
k,?
k
记估计的均方误差阵Pj?
k
^~
Xj=X-Xj?
kj?
~~tP=E（Xj,Xj）j?
k?
~最小均方误差估计就是指P为最小值，即Xj为X线形最小方差估计。
j?
kk离散线性kalman滤波的表达式
X=Z*X+G*Wkk-1k-1
Y=H*X+Vkkkk
其中,动态噪声{W}与测量噪声{V}是互不相关的0均值白噪声系列,即对kk
所有的k,j,系统噪声的统计特性为
EW=0,EV=0;
kk
W）=Qδ（k-j）kjkCOV（W
V）=Rδ（k-j）k,jkCOV（V
V）=0k,jCOV（W
1k,0,（k）,,,0k,0,
设系统初始状态的统计特性为
EX=μ00
tVar（X）=E[（X-μ）（X-μ）]00000且X与W，V都不相关,即0kk
COV（X,W）=COV（X,V）=00k0k
6
以上关于系统噪声的假设称为标准系统噪声假设。
初始值X,P,k=000
k=k+1
^^一步预报值X=GX?
kk-1k-1
tt=ZPZ+GQG预报误差P?
tt-1滤波增益K=PH（HPH+R）?
?
kkk-1kk-1
^^^最优滤波值X=ZX+K（Y-H*X）kk-1kkk
）?
k-1
滤波误差方差P=（I-KH）P?
kkkk-1
图四KALMAN的计算方法流程图
本题中的目标模型
7
X=X+Wkk-1k-1
Y=X+Vkkk
其中XY,WV是一维随机变量;
噪声{V}{W}与初始状态X满足标准系统假设;
kkk,kkK0
EX=μ00
）=P00Var（X附录二:
程序清单
数据处理程序:
%数据处理程序
ys1;
%ys1为去除斜三列后的数据ys1=fliplr（ys1）;
%左右转换
%修正算法程序
fori=1:
1:
24
forj=1:
13
ifi+j<
25
-j）.*ys1（i,j）;
szys（i,j）=（24-i）/（25-i%使用月数是j-1，不是j
else
szys（i,j）=0;
end
%修正千车故障数差分程序
dcf=szys;
%szys为修正后的数据
dcf1=zeros（24,13）;
fori=1:
21
if22>
i>
11
dcf1（i,1:
（23-i））=[dcf（i,1）,diff（dcf（i,1:
（23-i）））];
13）=[dcf（i,1）,diff（dcf（i,1:
13））];
ifdcf1（i,j）<
0
dcf1（i,j）=0;
%把0301批次数据去除
8
dcf2=dcf1;
%dcf1为差分后的数据dcf2（13,:
）=[];
dcf2;
%横向最小二乘法与纵向kalman滤波联合预测dcf3=dcf2;
%dcf2为去除0301批次后的数据fork=13:
23
%横向最小二乘法预测
y=dcf3（k,1:
（23-k））;
%多项式拟合
yy=conmin2（y）;
%纵向kalman滤波平滑
x=（1:
（k-1））;
y=[dcf3（x,（24-k））;
yy];
y=y'
;
%线性kalman滤波
zz（i）=conkalman1（y）;
dcf3（k,（24-k））=zz（i）;
end
dcf4=dcf3;
fori=12:
（24-i）
y=dcf4（（j+11+i-12）,1:
（14-j））;
x=1:
（j+i-2）;
y=[dcf4（x,（14-j））;
dcf4（（i+j-1）,（14-j））=zz（i）;
yy=fliplr（dcf4）;
%横向拟和与纵向Kalman滤波后的差分预测数据
functionxx=conkalman1（y）%线性kalman滤波程序
t=length（y）;
q=0.01;
r=0.01;
x0=0;
p0=100;
xx1=[];
k=[];
p=[];
9
t
ifi==1
k（i）=（p0+q）./（p0+q+r）;
xx1（i）=x0+k（i）.*（y（i）-x0）;
p（i）=r.*k（i）;
k（i）=（p（i-1）+q）./（p（i-1）+q+r）;
xx1（i）=xx1（i-1）+k（i）.*（y（i）-xx1（i-1））;
xx=xx1（end）;
functionyy=conmin2（y）%最小二乘曲线拟合预测
%输入为一行向量
%输出为预测的数
%数据的预处理
k=length（y）;
k;
%y=fliplr（y）;
%1-AGO的生成
forj=2:
y（j）=y（j-1）+y（j）;
ifk>
=1
%拟合实现的预测
p=polyfit（x,y,1）;
yy1=polyval（p,[x,x（end）+1]）;
%数据的后处理
%1-AGO
yy=yy1（k+1）-yy1（k）;
else
yy=0;
yy;
10
升级会员 