2011年全国大学生统计建模大赛论文一等奖第二名.pdf

I 不良产品返回数预测统计模型不良产品返回数预测统计模型 华南师范大学华南师范大学 莫祺、张先灯、杨鑫莫祺、张先灯、杨鑫 目录目录 摘 要 1 一、引言 2 二、数据来源及介绍 2 三、数据分析 3 四、模型一：联合多项分布概率模型 4 五、模型二：基于联合多项分布的 logistic 概率模型. 9 六、模型改进和解释. 14 七、两个模型的讨论与分析. 17 八、模型的不足与改进. 18 参考文献. 19 附件. 20 1 摘摘 要要 随着质量监管部门监管力度的加大和消费者维权意识的增强， 生产商必须更 好地履行质量三包的协议.然而产品的失效有一定的随机性，而生产商又想知道 该安排多少维修工人参与售后服务， 因此就有必要预测每一批生产产品的不良返 回数（失效产品的数目）.本文根据恩斯迈电子（深圳）有限公司数据部提供的 产品的一些月份的生产数据和相应的月份在一些月线的不良产品返回情况， 提出 了联合多项分布概率模型和基于联合多项分布的 logistic 概率模型两个模型，利 用 excel 和 sas 求出一段时间内的失效率，进而预测不良返回数.通过实例数据和 模型计算结果发现， 模型一较差， 通过加入生产总数和月份等协变量结合 logistic 回归得到模型二，预测结果有了很大的改进.另外，我们模型中的区间估计的区 间长度相对来说较短，验证的真实值落入的百分率不高，但只要做一点调整就可 达到超过 80%的准确率.也许是我们模型的拟合度不是特别高，也许是数据方面 的问题，有几个月线的不良返回数的偏差较大， 我们认为结合产品的专业知识来 寻找一些显著的协变量可能对模型会有更大的改进. 关键关键词词： 失效 不良返回数 多项分布 logistic 协变量 区间估计 2 一、一、引言引言 随着人们法律意识的增强，消费者的维权意识也越来越强.每年的 3.15 消费 者权益保护日报道出来的消费者投诉的事件也越来越多， 消费者对产品的质量也 越来越重视，国家也在逐渐完善质量监管体系.但是从概率的角度来说，谁也不 能保证某件产品的寿命.并且产品质检的时候不可能每个产品都来检测，一来耗 费人力物力，二来对有些产品是有毁坏性的.为了保证消费者的权益，商家就必 须对自己的产品承诺质量三包，尤其是电子产品，并且商家的售后服务对于商家 的声誉和销售都很重要 1. 每个厂家几乎都有售后服务部， 而售后服务需要投入多少维修人员合适也是 一个值得考虑的问题，换句话说就是要知道有多少产品有质量问题，需要返厂. 所以从生产商的角度来说， 预测一批产品有多少不良产品需要返厂维修就显得尤 为重要了 2. 二、二、数据来源及介绍数据来源及介绍 我们的数据来源于恩斯迈电子（深圳）有限公司数据部.恩斯迈电子（深圳） 有限公司为台湾股票上市公司微星科技股份有限公司（计算机主板、显卡、伺服 器全球前五大生产企业）投资于大陆的全资子公司，主要经营电子产品.而电子 产品的保修期一般都是三年， 因此有必要对一批产品出厂后三年内的不良产品的 返厂数进行统计和估计.考虑现实情况和经验，他们对一批产品的返厂数目一般 分了六个特殊时间段进行统计， 分别记为一月线不良返回数， 三月线不良返回数， 六月线不良返回数，十二月线不良返回数，二十四月线不良返回数，三十六月线 不良返回数.而其实一月线的返回数是指这批产品出厂后的三个月内的不良返回 数， 三月线是指这批产品出厂后的五个月内的不良返回数，六月线不良返回数指 的是八个月内的返回数，同理十二月线、二十四月线、三十六月线都要往后推迟 两个月，这样的话就基本上可以统计保修期内返回数，反过来又可以判断这批产 品的质量情况.他们把 2007 年一月至 2010 年十一月的部分返回数据统计出来了. 由于最初拿到数据在 2011 年二月统计完的， 所以 2010 年十月和十一月我们只能 统计这两批产品的一月线的返回数，2010 年六、七、八月可以得到一月线和三 月线的返回数，2010 年一月至六月可以得到一月线、三月线和六月线的返回数， 2009 年一月至十二月可以得到一、三、六、十二月线的返回数，2008 整年只是 不能得到三十六月线的不良返回数，而 2007 年就可以得到各个月线的返回数据. 而作为生产商来说，就想知道 08 年各月的三十六月线的返回情况，09 年的二十 四月线及以后的返回情况，10 年的一月至六月的十二月线及以后的返回情况， 等等. 后来，我们还得到该公司 MB、VGA、DT-BB、EPS 四种电子产品的各月三 十六个月线的不良数的返回情况，这对我们的预测提供了很大的帮助. 3 三、数据分析三、数据分析 这里的数据预测其实就是要求预测每个月相应的各个月线的不良数的返回 情况， 可以看做一个纵向数据的类型，每个月的各月线数据组合在一起就变成了 一个联合纵向数据类型.尽管纵向数据已经有很多研究方法，尤其在心理学的研 究中占据着重要的地位， 比如时间序列的方法， 混合效应模型， 潜变量增长模型， 多层线性模型等等 3.对于时间序列来说，需要的时间数据多些效果才会比较好， 一般来说至少需要 50 个，而我们这里只提供 6 个时间点的数据，显然不合适. 后面几种方法都是基于回归的思想发展起来的方法，本文根据数据的特点，将联 合多项分布与 logistic 回归结合起来，来预测产品在每个时间区间的失效率，进 而来估计产品的不良返回数.实例的数据结构可以简化为下表的形式： 一类月 二类月 三类月 四类月 五类月 六类月 生产数 1 n 2 n 3 n 4 n 5 n 6 n 一月线返回数 11 n 21 n 31 n 41 n 51 n 61 n 三月线返回数 12 n 22 n 32 n 42 n 52 n 六月线返回数 13 n 23 n 33 n 43 n 十二月线 返回数 14 n 24 n 34 n 二十四月线 返回数 15 n 25 n 三十六月线 返回数 16 n 表 1 表 1 中一类月的数据类型就和 07 年的各月的数据类型一样的，知道各个月 线的不良返回数情况，二类月与 08 年各月的数据类型一样的，仅仅不知道三十 六月线的返回情况，同理知，三类月与 09 月的一样，四类月与 10 年一月至六月 一样，五类月与 10 年七、八、九月一样，六类月与 10 年十月和十一月一样.表 中生产数表示这个月生产产品的数目， 各个月线的返回数与数据来源中的介绍是 一样的.表中 i n表示第 i 类月生产的产品数， ij n表示第 i 类月生产的产品数在 j 月 线的不良产品返回数( ,1,2,6)i j L，而作为厂家来说，就希望将这个表中那些 空白的部分预测出来. 4 四、模型一四、模型一：联合多项分布概率模型联合多项分布概率模型 1、建立模型、建立模型和模型和模型求解求解 假设所有产品的生产工艺和技术都是一样的， 对于一个产品是否失效可以看 作一个随机事件，根据题中数据的分布情况，我们可以假设每件产品在出厂后三 个月内月损坏的概率为 1 p，第四个月到和五个月之间损坏的概率为 2 p，第六个 月至第八个月损坏的概率为 3 p，第九个月至第十四个月之间损坏的概率为 4 p， 第十五个月至第二十六个月之间损坏的概率为 5 p，第二十七个月至第三十八个 月之间损坏的概率为 6 p.那么对于一类月的数据的分布情况，每件产品在一月线 时间内损坏的概率为 1 p，第四个月到和五个月之间损坏的概率为 2 p（一月线内 没有损坏， 三月线内损坏的概率为 2 p） ， 第六个月至第八个月损坏的概率为 3 p（三 月线内没有损坏，六月线内损坏的概率为 3 p） ，第九个月至第十四个月之间损坏 的概率为 4 p（六月线内没有损坏，十二月线内损坏的概率为 4 p） ，第十五个月至 第二十六个月之间损坏的概率为 5 p（十二月线内没有损坏，二十四月线内损坏 的概率为 5 p） ，第二十七个月至第三十八个月之间损坏的概率为 6 p（二十四月线 内没有损坏，三十六月线内损坏的概率为 6 p） ，三十六月线以内不坏的概率为 123456 1pppppp.同理我们也可知道二类月、 三类月、 四类月、 五类月、 六类月的概率分布情况.我们可以把表 1 换成 一类月 二类月 三类月 四类月 五类月 六类月 生产数 1 n 2 n 3 n 4 n 5 n 6 n 一月线返回数 11 x 21 x 31 x 41 x 51 x 61 x 三月线增加的 返回数 12 x 22 x 32 x 42 x 52 x 六月线增加的 返回数 13 x 23 x 33 x 43 x 十二月线增加的 返回数 14 x 24 x 34 x 5 二十四月线增加 的返回数 15 x 25 x 三十六月线增加 的返回数 16 x 表 2 其中 11ii xn ， ,1( 2) ijiji j xnnj ，这样做之后相当于将原始数据中的累 积不良返回数换成了各个时间段增加的不良返回数.这样的话，我们可以把一类 月的生产数与各月线增加的返回的情况看做一个总数为 1 n七项分布， 并且一类月 的数据的似然函数 1315161111213141516111214 11123456123456 X(1) xxxnxxxxxxxxx LC pppppppppppp 1 （ ） 其中 1 1 1112131415161111213141516 ! !()! n C xxxxxxnxxxxxx 111213141516 X(,)xxxxxx 1 同理我们知道二类月的数据可以看做一个六项分布，相应的似然函数为 232522122232425212224 2221234512345 X(1) xxnxxxxxxxx LC pppppppppp （） 而 2 2 212223242522122232425 ! !()! n C xxxxxnxxxxx 2122232425 X(,)xxxxx 2 类似地一直到六类月，分别是五项分布，四项分布，三项分布，二项分布，我们 也可以得到相应的似然函数 3456 ,L L L L 4.利用上表中所有数据建立的似然函数 为 65432 12345 1111116 1111213141516 22122232425 331323334 123456123456 123456 12345 1234 (,) (1) (1) (1) (1 iiiii iiiii xxxxx x nxxxxxx nxxxxx nxxxx L XXXXXXCpppppp pppppp ppppp pppp 4414243 55152661 123 121 ) (1)(1) nxxx nxxnx ppp ppp 这里面有 123456 ,p ppppp六个参数，对上述似然函数取对数，然后分别对 123456 ,p ppppp求偏导可以得到的 123456 ,p ppppp的极大似然估计的递推 关系式 6 (1)(2)(3)(4)(5)(6) 112131415161112131415161 1 123456123456 xxxxxxn pn pn pn pn pn p p nnnnnnnnnnnn 1222324252 21 111221331441551 (1)(2)(3)(4)(5) 1222324252 1 (1)(2)(3)(4)(5) 1121314151 (1) ()()()()() (1) 11111 xxxxx pp nxnxnxnxnx n pn pn pn pn p p npnpnpnpnp （） （） （） （） （） 13233343 312 11112221223313244142 (1)(2)(3)(4) 13233343 12 (1)(1)(2)(2)(3)(3)(4)(4) 112212312412 (1) ()()()() (1) 1111 xxxx ppp nxxnxxnxxnxx n pn pn pn p pp nppnppnppnpp （-） （-） （-） （-） 142434 4123 111121322122233313233 (1)(2)(3) 142434 123 (1)(1)(1)(2)(2)(2)(3)(3)(3) 112321233123 (1) ()()() (1) 111 xxx pppp nxxxnxxxnxxx n pn pn p ppp npppnpppnppp （-） （-） （-） 1525 51234 111121314221222324 (1)(2) 1525 1234 (1)(1)(1)(1)(2)(2)(2)(2) 1123421234 (1) ()() (1) 11 xx ppppp nxxxxnxxxx n pn p pppp nppppnpppp （-） （-） 16 612345 11112131415 (1) ( ) 16 12345 (1)(1)(1)(1)(1) 112345 (1) () (1) 1 iji j i x pppppp nxxxxx x n p pppppp npppppn 其中 （-） ( ) i j p为单独用第 i 类月的数据对 j p的似然估计.当(1,2,6) i n i L很大的时候， 并且(1,2,6) j pj L很小的时候， 123456 ,p ppppp的估计值可以近似简化为 ( 1 )( 2 )( 3 )( 4 )( 5 ) 1222324252 112 12345 , n pn pn pn pn p ppp nnnnn (1)(2)(3)(4)(1)(2)(3) 13233343142434 34 1234123 , n pn pn pn pn pn pn p pp nnnnnnn 7 (1)(2)(1) (1) 152516 566 121 , n pn pn p ppp nnn 此时 j p的估计 j p 可以看做是各 个月估计 ( ) i j p的加权平均值，相比较前面递推关系式，可以简便很多计算，而 且估计也不会差很远.并且由于 i n一般都很大， j p比较小，由多项分布的性质知， (,)( ,1,2,6) mn r ppm nmn L且几乎趋于 0,也就是说 m p 与 n p 可以看做是不相关 的，也即是相互独立的.一旦把 123456 ,p ppppp估计出来，我们就可估计表中 空着的数据 ijij xn p ， 比如二类月的三十六月线增加的返回数 2626 xn p ，因 此二月的三十六月线的返回数的估计就为 2526 xx . 2、区间估计、区间估计 我们知道点估计对样本的依赖性很大，而且有很大的误差，并且没办法衡量 偏差程度.区间估计则可以按一定的可靠性对待估参数给出一个区间范围，因此 我们有必要估出不良返回数的置信区间. 多项分布的边际分布是二项分布，所以(),()(1) ijijijijj E xn p D xn pp 当 i n和 ij x的个数比较多时，我们可以近似认为(1) ijijijj xn pn ppN(,同理我们 也知道 (1 ) jj jj i pp pp n ） N(,，因此(1) ijijijj n pn pn pp N(,，所以对于概率 ()()()() 22 ijijijjjijijjjij cc P xn pcP xnpnpn pcP xnpP npn p 00 0 ()0.005()0.005 22 ijjjij cc cP xnpP npn p 如果我们找到某个 使得， 所以 0 ()0.01 ijj P xn pc ，因此 ij x的 99%置信区间为 00ijij n pcn pc ， 查表知 0 5.62(1) ijj cn pp .，但是为了方便我们就取 0 6(1) ijj cn pp 来得 到置信度更高的区间.例如 26 x的置信度为 99%的置信区间为： 2626626266 6(1)6(1)n pn ppn pn pp ， 因此二类月中三十六月线不良返回数的 99%置信区间为 8 2 5262662 526266 6(1)6(1)xnpnppxnpnpp ， 3、模型结果、模型结果 模型一中我们并没有假定 123456 ,p ppppp分布情况， 对 123456 , , , , ,p p p p p p 的估计完全依赖于样本，借助于 47 个月的部分数据情况建立的似然函数 1247 (,.)L X XX取对数，再 123456 ,p ppppp分别求偏导数，借助于 excel 里面 的菜单操作计算功能可以得到估计值 123 0.002846,0.003875,0.005573ppp , 4 0.009478p , 56 0.01655,0.008133pp .而后来我们得到了三个月之后的一些 不良返回数据，其实就是 2010 年十月和十一月的三月线，2010 年七月至九月的 六月线， 2010 年一月至三月的十二月线， 2009 年一月至三月的二十四月线， 2008 年一月至三月的三十六月线，我们刚好用这些数据来验证的我们的模型，得到如 下表的结果 月份及月线 估计值 置信区间 真实值及偏差 08 年一月三十六月线 53356 （52045，54691） 60407，-7039 08 年二月三十六月线 49346 （48074，50617） 55072，-5726 08 年三月三十六月线 77703 （761181，79226） 91833，-14130 09 年一月二十四月线 34909 （33826，35992） 38119，-3210 09 年二月二十四月线 42501 （41245,43757） 36450, 6051 09 年三月二十四月线 52178 （50791，53565） 44800， 7378 10 年一月十二月线 25066 （24131，26001） 24142， 924 10 年二月十二月线 21012 （20148，21876） 20172， 840 10 年三月十二月线 25120 （24256，25984） 24550， 570 10 年七月六月线 6291 （5782,6800） 5323， 968 10 年八月六月线 10361 （9710，11012） 9905， 456 9 10 年九月六月线 9351 （8711,9991） 8316， 1035 10 年十月三月线 3761 （3422,4100） 2672， 1089 10 年十一月三月线 6104 （5683，6525） 6402， -298 表 3 注：表中的置信区间是 99%置信区间，偏差=估计值真实值.有阴影月份的月线表示真实值 落在我们预测的置信区间内. 从表 3 的验证结果可以看出，我们的预测效果并不怎么好，14 个月份的月 线预测结果只有 5 次落在置信区间内， 准确率不到 40%， 并且不在置信区间内的 月份的偏差大多数都比较大，结果比较糟糕.但是我们观察原始数据中的累积不 良率 AFR （ ( ) t AFR t 时刻累积的不良数 生产总数 ）发现， 各个月在同一个时间点的 AFR 差异还是很大的，譬如一月线的 AFR,07 年 2 月的一月线 AFR 才 0.12%,而 08 年 5 月的一月线 AFR 却为 0.53%，可见 AFR 在一月线这个时间点的时间点的跨度 是很大的，再看三月线的 AFR，最小的才 0.35%，最大的达到 1.1%，后面的月 线也出现类似的结果，而我们的模型中在产品的失效率只与时间段或点有关，其 实也就是只与月线有关，只要是同一月线我们认为失效率都是一样的.而样本数 据告诉我们，同一月线在不同的月份的失效率差异很大，所以我们的模型单用纵 向时间来来估计失效率是不好的， 因此我们还要寻找其他的因素来估计在某个时 间段的失效率，因而我们就提出了模型二. 五五、模型、模型二：二：基于基于联合多项分布联合多项分布的的 logistic 概率模型概率模型 1、模型建立、模型建立 由模型一知，我们数据中失效率不仅仅与纵向时间有关.虽然从可靠性理论 来说 5，一件产品的失效主要是与使用时间有关，但是也不能排除某些协变量会 对我们的失效率有影响.由于在某一个月份时对 123456 ,p ppppp的估计我们是 用对应的不良返回数除以该月的生产总数，所以我们就有必要考虑生产总数对 123456 ,p ppppp的影响，所以就把生产总数当作协变量进我们的模型中.由于 123456 ,p ppppp的范围都在 01，而生产总数的值与他们之间的相差很大，因 此我们考虑 logistic 模型，对 123456 ,p ppppp作 logistic 变换，从而可以使得它 们的范围变成负无穷到正无穷 6. 10 对于我们要研究的数据，有 K=47 个月的生产的产品要来估计不良返回数. 假定每个月的观察的时间 01 0(6) m tttm.由给定的六个月线的时 间关系知， 123456 3,5,8,14,26,38,tttttt各个月累积的不良返 回数和新增的不良数可以表示为(,):1,6;1,47 kjkj dxjkLL，其中 1k 表示 07 年一月，2k 表示 07 年二月，以次类推，47k 表示 10 年十 一月， kj d表示第 k 月在 j t时刻累积的不良返回数， kj x表示第 k 月生产的产品在 时间段 1 (, ) ijj ttI 期间即将失效的产品个数，用 kj p表示失效率，可以解释为 第 k 月生产的一件产品在时间段 j I的失效率，同时我们用 k n表示第 k 月生产的 产品总数，那么第 k 月在 j t时刻还没有失效的产品总数可以表示为 k jkk j rnd， 根据模型一中我们对数据的纵向分析易知， 在时间段 j I期间第 k 月中将要失效的产品服从二项分布，即 ,1 (;) kjk jkj xB rp . 我们建立 logit 链接函数 log()log()() 1 kj kjkjkj kj p it pmg te p 7 (1) 这里g是具有未知参数的一个已知函数， k m是模型的协变量，可能的向量取 值 12 (1,) kkk mmm， 012 (,) 是对于于协变量的回归系数向量.而 kj e是时间区间 j I内第 k 月产品失效率的随机效应.设 1 (,) , kkkm eeeL其中 1,2.kKL 假定这些随机效应 kj e在同一个月中相关， 而在任两个不同的月中 相互独立.设 12 ( ,) m Adiag a aaL， 这里 1/2 1 ()(1,) jjj attjm L是 连续的时间长度的平方根.我们用 m 维多元正态分布来建立每个月的相关结构， 满足 2 (0,) km eN 这里是m m矩阵，定义为 11 21 2 123 1 1 1 m m mmm AA L L MMMM L ( )gt 是未知参数确定的时间函数，( )gt 的选定依赖于这类失效数据的反应 时间模型.最简单的形式是线性，即 1 ( )g tt ，没有常数项的是为了保证模 型 的 可 识 别 性 . 对( )gt 的 选 择 还 有 很 多 其 它 的 形 式 . 如 二 次 式 2 12 ( )g ttt ，甚至更复杂的.但是( )gt 并不见得越复杂越好，比如 Waterloo 大学生物系进行的一系列的大种系动物的毒理实验中多元数据分析来 看 ， 他 们 用 贝 叶 斯 的 方 法 来 选 择( )gt 时 发 现 ， 最 佳 的 模 型 是 12 ( )g tt 10 11 . 对于此时模型的区间估计，我们的思路和模型一的原理是一样的.这里也从 多项分布和二项分布之间的关系出发，结合正态分布来寻找的置信区间.由于不 同月份生产的产品在同一个时间段内的失效率是不同的， 所以模型一区间估计中 的 j p要换成(1,2,) ij p iKL. 2、模型结果、模型结果 我们同样利用模型一的数据来拟合模型，并且也同样用那 14 个返回数据来 验证模型， 由于此时的数据不多， 我们只加入各月的生产总数作为协变量，( )gt 我们就选择了关于时间的一次函数，利用 sas 编程 6 8得到的结果如下： 从方差分析的 p 值小于 0.0001 知，模型拟合数据时显著性有效的，而对于 未知参数估计的 p 值也都显著性不为 0，所以我们的模型为 12 log()3.77912770.364566.70140 1 kj kj kj p Ent p (2) 也即是 ( 3.77912770.364566.70140) 1 1 kj kj Ent p e (3) 这里 kj p表示第 k 月产品的在 j I时间段的失效率， k n表示第 k 月的生产的 产品总数.而 j t代表的月线，也就是分布为一月线、三月线、六月线、十二月线、 二十四月线、三十六月线，但是在我们的模型中的取值分别为 1,、2、3、4、5、 6.我们也将 j t分别取 3、5、8、14、26、38 时的取值来拟合模型，但是最后的预 测效果没有我们的模型好. 为了直观地评估模型的拟合度，模型的学生化残差的散点图如下： y = -6.7014 +378E-9n +0.3646 t N 206 Rsq 0.7348 AdjRsq 0.7322 RMSE 0.3376 Studentized Residual -4 -3 -2 -1 0 1 2 3 Predicted Value -6.25-6.00-5.75-5.50-5.25-5.00-4.75-4.50-4.25-4.00 图 1 从图 1 我们可以看出大多数标准化残差点都在两倍标准差之内， 只有少数点 落在之外，此时对于十四个验证的数据的估计如下： 月份及月线 估计值 置信区间 真实值及偏差 08 年一月三十六月线 62323 （60868，63778） 60407，1916 08 年二月三十六月线 57096 （55702，58489） 55072，2024 13 08 年三月三十六月线 92920 （91151，94690） 91833，1087 09 年一月二十四月线 31554 （30539，32568） 38119，-6565 09 年二月二十四月线 36815 （35682,37948） 36450, 365 09 年三月二十四月线 47002 （45723，48282） 44800， 2202 10 年一月十二月线 23437 （22528，24346） 24142， -705 10 年二月十二月线 19143 （18321，19965） 20172， -1029 10 年三月十二月线 23529 （22618，24440） 24550，-1021 10 年七月六月线 5689 （5238,6139） 5323， 366 10 年八月六月线 10051 （9453，10650） 9905， 146 10 年九月六月线 8992 （8426,9558） 8316， 676 10 年十月三月线 3355 （3009,3702） 2672， 683 10 年十一月三月线 5956 （5496，6417） 6402， -446 表 4 从表 4 知， 有 7 个月的月线的真实值在预测的区间内， 但仔细观察你会发现， 10 年一月、二月的十二月线的的真实值与预测的置信的区间的置信上限只不过 相差 100 左右，只要我们稍作调整就可以落入置信区间内.而对于 08 年一月、二 月的三十六月线与 09 年三月二十月线，它们预测值与真实值之间的偏差大约都 在 2000 左右，这个偏差还不到真实值的 1/20，在要求不是很严的前提下还是可 以接受的.而对于 09 年一月二十四月线和 10 年十月的三月线感觉就不是很让人 接受， 但仔细观察原始数据发现， 09 年一月二十四月线的不良返回数为 20274,09 年二月二十四月线的不良返回数为 22820，而 09 年一月三十六月线的真实不良 返回数为 38119,09 年二月三十六月线的真实不良返回数为 36450， 二十四月线是 09 年一月比二月时的返回数少了 2546，但是在三十六月线时一月的返回数比二 月的返回数却多了 1669，所以 09 年一月的数据变异性较大.但总体来说，比模型 一的结果还是好了很多. 14 六、模型六、模型改进和解释改进和解释 1、模型改进、模型改进 由于后来我们得到了这个产品（MB）一些月份三十六个月线的返回数据， 所以此时我们就用了每个月一月线、二月线、一直到三十六月的数据来拟合模型 二.由于这时的数据多了很多，所以此时的协变量只取每个月的生产总数显然不 合适.我们就想是否具有季节性，因此在协变量中加入了月份，但是月份的影响 是周期性的，故构造了周期为 12 的三角函数对失效率的影响，模型为： log ()log()( ,) 1 kj kjkjkkj kj p it pmg t ne p （4） 此时的(1,cos(),sin() 66 kk mnss 其中mod( ,12)sk， 另外我们还加 入了生产总数与月线的交互项 kj n t，拟合数据得到的结果如下： 从 sas 运行的结果可以看出，模型的检验和系数的检验都可以通过.其中 sn 表示 的是sin() 6 s ，cn 表示的是cos() 6 s ，故此时的模型可写为： log()1.44029170.059362.5182380.04607sin()0.18730cos()6.33065 166 kj kjkj kj p EntEntss p 也即是 ( 1.44029170.059362.5182380.04607sin() 0.18730cos() 6.33065) 66 1 (5) 1 kjkj kj EntEntss p e 这里 kj p表示第 k 月产品的在 j I时间段的失效率， k n表示第 k 月的生产的 15 产品总数.而 j t代表的月线，也就是分别为一月线、二月线、三月线、L、三十 六月线，但是在我们的模型中的取值分别为 1,、2、3、L、36.而 s 表示的是将 k 对 12 取余，即mod( ,12)sk，得到的标准化残差图： y = -6.3307 +144E-9n -0.0594 t +25E-9 nt +0.0461 sn +0.1873 cn N 1170 Rsq 0.5108 AdjRsq 0.5087 RMSE 0.3645 Studentized Residual -8 -6 -4 -2 0 2 4 Predicted Value -8.00-7.75-7.50-7.25-7.00-6.75-6.50-6.25-6.00-5.75 图 2 从标准化的残差图 2 可以看出，数据几乎都落在 2 倍标准差范围之内，只有 极少落在之外.我们后来也做了将标准差之外点剔除再来拟合模型，但是由于偏 差数据太少的缘故，对我们结果的影响很小. 月份及月线 估计值 置信区间 真实值及偏差 08 年一月三十六月线 60720 （60541，60900） 60407，313 08 年二月三十六月线 55229 （55069，55389） 55072，157 08 年三月三十六月线 92128 （91893，92362） 91833，295 09 年一月二十四月线 38171 （37993，38350） 38119，52 09 年二月二十四月线 36873 （36646,37100） 36450, 423 09 年三月二十四月线 45436 （45169，45702） 44800， 636 10 年一月十二月线 24738 （24469，25006） 24142， 696 16 10 年二月十二月线 20189 （19957，20420） 20172， 17 10 年三月十二月线 24637 （24416，24858） 24550，87 10 年七月六月线 5411 （5250,5572） 5323， 88 10 年八月六月线 9489 （9264，9714） 9905， -416 10 年九月六月线 8316 （8297,8759） 8316， 212 10 年十月三月线 3156 （2944,3369） 2672， 484 10 年十一月三月线 6446 （6170，6722） 6402，44 表 5 从表 5 可以看出，此时由于我们月份之间的跨度比较小，所以此时置信区间 的长度也比较短，所以此时的偏差总体上也比之前小了很多.08 年一月、三月的 三十六月线和 09 年二月二十四月线的真实值与置信区间的界限最近也就相差一 百多，经过小小调整就可以了.从上面的偏差知，偏差最大的也不过六百多，离 落入置信区间最大也不过三百左右，如果对精度不是要求很高的话，完全可以通 过增大区间的长度来落入置信区间内. 2、模型解释模型解释 对于上面式（5）中关于 kj p的关系式，当 k n 变大时，由于 k n 前面的系数 是负的， kj p也将会变大.这很容易理解，因为当一个月的生产总数增加了，而此 时工人的生产任务就加大了，而工人的精力有限，从而生产效率降低，因而就会 有更多的次品.当 j t 变大时，由于 j t 前面的系数是正的， kj p反而会减小.从可靠 性理论的角度来说，一个产品的失效的概率分布通常可以看做一个指数函数，由 于指数函数的分布函数( )1 t F te ， 我们可假设一个相隔一个时间间隔的概 率函数 (1) ( )( )(1) tt Q tF tF tee ， 而求导得 (1) ( )() tt Q tee 所以( )Q t是一个关于 t 的减函数，也就是说随着时间的增大，相等的时间间隔 内失效的概率是减小的，这刚好也我们的模型是吻合的.而对于月份的影响也是 很容易解释的，由于每个月之间的气候、湿度等等自然因素的不同，对产品的保 存也会造成一定的影响，故月份也是有影响的，只是说相对来说小点而已. 17 七七、两个模型的讨论与分析、两个模型的讨论与分析 两个模型从预测的方式来说可以说方式是一样的， 都是先预测一段时间内的 失效率， 再用失效率乘以生产总数得到增加的不良返回数，但是预测失效率的方 法是不一样的，并且引入的参数也不一样.模型一中我们的模型假设是在每一件 生产产品在出厂后的相同的一段时间内的概率是相同的，只依赖产品的出厂时 间，从而借助多项分布来构造样本的似然函数来估计各个区间的失效率.但是例 子验证的结果并不理想， 区间准确率约为 36%， 而其他落在区间外的月份的月线 的预测都离区间比较远，所以效果并不好.而模型二就引入了协变量来估计失效 率.除了考虑出厂时间，我们还考虑了生产总数的影响和月份的季节性影响，甚 至还考虑他们的之间的交互， 准确率大约维持在 50%多， 但是我们发现还有几个 月的真实值离我们的预测区间只有 100 左右的偏差， 我们的预测区间本来是就比 较短，只要区间稍微放大，我们就可以使得准确率超过 80%.对于改进的模型来 说，偏差最大也就 600 左右，也就是说如果我们把置信区间的长度放大到 1200 多，那我们就有可能全部预测对.此时的预测都是基于上一月线已知的不良数加 上新增的不良数来预测下一月线的不良数，因此我们就想能不能跳跃性地预测， 预测之后几个月线的.因此我们用改进的模型基于二十四月线的数据来直接预测 三十六月线的数据，也就是我们要求二十四月线至三十六月线之间增加的不良 数， 我们利用模型分别算出了二十四月线至二十五月线的增加数，二十五月线至 二十六月线的增加数，L，三十五月线至三十六月线的增加数，再加这十二个 间隔的增加数加起来预测二十四至三十六月线的增加数， 我们计算出来的结果发 现全部偏差很大，甚至还比不上我们之前就只用一、三、六、十二、二十四、三 十六月线的数据拟合的模型， 也就是说我们这个模型只适合预测一个时间点的数 据， 并且这个数据只与上一个时间点数据有关的数据类型，与更前面的的数据关 联不大，因此也不适合作跳跃性的时间预测. 对于模型二中存在小部分数据的学生化残差不在两倍标准差以内的数据， 我 们剔除了这些点再来拟合模型，结果发现对最后的结果影响不大，也就那些数据 对模型的贡献不怎么大.对于模型中出现的生产总数与月线的交叉项，我们也用 过不含交叉项的模型去拟合数据，结果和我们上面的差不多，都只是有小小的变 动，所以也可以考虑用不含交叉项的模型拟合. 从实例的结果来看，尽管模型二的结果相对于模型一的结果来说有所改进， 但是结果还不是特别满意，也只能说还可以，虽然改进的模型中预测的最大偏差 不过 600 多，但是我们仔细观察数据，发现从上一月线到下一月线增加的不良返 回数也不过就是 2000 左右， 这样来说 600 也还是比较大的偏差.对于此例的偏差， 我们分析可能有以下两方面的原因： （1）一方面可能来自于我们的模型，可能是 我们的模型不是太适合这个数据，可能是有某些重要的协变量我们没有找到，这 可能涉及到电子产品制造的某些专业知识，也正是我们所欠缺的.（2）另一方面 可能来自于数据，我们的数据来源于回收的不良产品数， 而我们的模型是应用于 产品的失效前提下，这两者是不等同的.因为有些产品失效了，但是可能却没有 及时回收.有些人怕麻烦或者比较忙没有及时将自己的失效的产品返厂，可能就 不返厂了或者拖到下个月、下下个月甚至更晚返厂.这样的话就会导致数据会有 偏差. 18 八、八、模型的不足与改进模型的不足与改进 本文结合联合多项分布和 logistic 回归提供了一种解决多个联合不完全纵向 数据的新的方法，虽然从实例来说不是特别理想，但是也不失为这种数据类型可 以选择的一种分析方法.我们是先通过预测失效率，再来预测产品失效的个数， 由于失效率比较小，相应的标准差也就小，所以我们求得的置信区间的长度不是 很大，也就说更强调精度，应用价值更好.我们两个模型的区间估计几乎是用一 样的方法，模型一中的模型假设某一段时间的实效率是与生产总数是无关的，而 模型二中失效率是与生产总数是有关系的，而考虑失效率一般都很小，即使有生 产总数的影响，偏差也不会大得很多，本来我们的区间就不长，所以我们认为是 可以接受的，当然从假设模型的分布得到模型二的区间估计是更有说服性. 我们的模型二中是借助于 logistic 模型来估计失效率，也许我们还可以借助 于多项式曲线回归、Gompertz 曲线回归、二项型指数曲线回归、三项型指数曲 线回归等等更复杂的模型来拟合失效率，或许可以得到一些更好的模型.另外， 从改进的模型的解释部分可知，我们的模型一定程度可以用指数分布来解释，然 而我们知道在产品寿命研究还有一个更复杂点的威布尔分布， 如果我们能构造一 个模型最终可以用威布尔分布来解释，可能会有意想不到的结果 6 9.对于协变 量的选取，我们模型选了生产总数和月份，可能还有某些重要的协变量我们没有 找到， 结合电子制造和工程学的相关知识对我们寻找协变量可能也会有很多的帮 助，由于得到数据的时间相对来说比较短，我们没有时间去学习这方面的知识， 但应该是一个可以考虑的方面. 19 参考文献参考文献 1张洁，产品质量责任的归责原则，中国商界，2010(4). 2陈琳、柳鹏飞，产品维修质量的综合评估和计算方法，质量与可靠性，2005 （1）. 3刘红云、孟庆茂，纵向数据分析方法，心里科学进展，2003，11（5）. 4茆诗松、程依明、濮晓龙，概率论与数理统计教程，高等教育出版社，2007. 5曹晋华、程侃，可靠性数学引论，.高等教育出版社，2006. 6胡良平，sas 统计分析教程，电子工业出版社，2010. 7方积乾、陆盈，现代医学统计学，人民卫生出版社，2002. 8朱世武，sas 编程技术教程，清华大学出版社，2007. 9薛仁言，威布尔模型族，科学出版社，1998. 10Chen, M.-H., Dey, D.K.,and Sinha,D.(2000).Bayesian analysis of multivariate mortality data with large families. Applied Statistics, 49, 129-144. 11OHara Hines, R.J.(1989). Some methods for the analysis of texicological mortaility data grouped over time. Unpublished PHD Thesis, Department of Statistics and Actuarial Science, University of Waterloo, Canada. . 20 附件附件 程序一 data chengxu1; array a7,47; array b7,47; do i=1 to 7; do j=1 to 47; input x; a(i,j)=x; end; end; do i=1 to 47; b(1,i)=a(1,i); b(2,i)=a(2,