数学建模论文-不良产品返回数预测统计模型.doc

不良产品返回数预测统计模型华南师范大学 目录摘 要1一、引言2二、数据来源及介绍2三、数据分析3四、模型一：联合多项分布概率模型4五、模型二：基于联合多项分布的logistic概率模型9六、模型改进和解释14七、两个模型的讨论与分析17八、模型的不足与改进18参考文献19附件2027摘 要随着质量监管部门监管力度的加大和消费者维权意识的增强，生产商必须更好地履行质量三包的协议.然而产品的失效有一定的随机性，而生产商又想知道该安排多少维修工人参与售后服务，因此就有必要预测每一批生产产品的不良返回数（失效产品的数目）.本文根据恩斯迈电子（深圳）有限公司数据部提供的产品的一些月份的生产数据和相应的月份在一些月线的不良产品返回情况，提出了联合多项分布概率模型和基于联合多项分布的logistic概率模型两个模型，利用excel和sas求出一段时间内的失效率，进而预测不良返回数.通过实例数据和模型计算结果发现，模型一较差，通过加入生产总数和月份等协变量结合logistic回归得到模型二，预测结果有了很大的改进.另外，我们模型中的区间估计的区间长度相对来说较短，验证的真实值落入的百分率不高，但只要做一点调整就可达到超过80%的准确率.也许是我们模型的拟合度不是特别高，也许是数据方面的问题，有几个月线的不良返回数的偏差较大，我们认为结合产品的专业知识来寻找一些显著的协变量可能对模型会有更大的改进.关键词： 失效 不良返回数 多项分布 logistic 协变量 区间估计 一、引言随着人们法律意识的增强，消费者的维权意识也越来越强.每年的3.15消费者权益保护日报道出来的消费者投诉的事件也越来越多，消费者对产品的质量也越来越重视，国家也在逐渐完善质量监管体系.但是从概率的角度来说，谁也不能保证某件产品的寿命.并且产品质检的时候不可能每个产品都来检测，一来耗费人力物力，二来对有些产品是有毁坏性的.为了保证消费者的权益，商家就必须对自己的产品承诺质量三包，尤其是电子产品，并且商家的售后服务对于商家的声誉和销售都很重要.每个厂家几乎都有售后服务部，而售后服务需要投入多少维修人员合适也是一个值得考虑的问题，换句话说就是要知道有多少产品有质量问题，需要返厂.所以从生产商的角度来说，预测一批产品有多少不良产品需要返厂维修就显得尤为重要了.二、数据来源及介绍我们的数据来源于恩斯迈电子（深圳）有限公司数据部.恩斯迈电子（深圳）有限公司为台湾股票上市公司微星科技股份有限公司（计算机主板、显卡、伺服器全球前五大生产企业）投资于大陆的全资子公司，主要经营电子产品.而电子产品的保修期一般都是三年，因此有必要对一批产品出厂后三年内的不良产品的返厂数进行统计和估计.考虑现实情况和经验，他们对一批产品的返厂数目一般分了六个特殊时间段进行统计，分别记为一月线不良返回数，三月线不良返回数，六月线不良返回数，十二月线不良返回数，二十四月线不良返回数，三十六月线不良返回数.而其实一月线的返回数是指这批产品出厂后的三个月内的不良返回数，三月线是指这批产品出厂后的五个月内的不良返回数，六月线不良返回数指的是八个月内的返回数，同理十二月线、二十四月线、三十六月线都要往后推迟两个月，这样的话就基本上可以统计保修期内返回数，反过来又可以判断这批产品的质量情况.他们把2007年一月至2010年十一月的部分返回数据统计出来了.由于最初拿到数据在2011年二月统计完的，所以2010年十月和十一月我们只能统计这两批产品的一月线的返回数，2010年六、七、八月可以得到一月线和三月线的返回数，2010年一月至六月可以得到一月线、三月线和六月线的返回数，2009年一月至十二月可以得到一、三、六、十二月线的返回数，2008整年只是不能得到三十六月线的不良返回数，而2007年就可以得到各个月线的返回数据.而作为生产商来说，就想知道08年各月的三十六月线的返回情况，09年的二十四月线及以后的返回情况，10年的一月至六月的十二月线及以后的返回情况，等等.后来，我们还得到该公司MB、VGA、DT-BB、EPS四种电子产品的各月三十六个月线的不良数的返回情况，这对我们的预测提供了很大的帮助.三、数据分析这里的数据预测其实就是要求预测每个月相应的各个月线的不良数的返回情况，可以看做一个纵向数据的类型，每个月的各月线数据组合在一起就变成了一个联合纵向数据类型.尽管纵向数据已经有很多研究方法，尤其在心理学的研究中占据着重要的地位，比如时间序列的方法，混合效应模型，潜变量增长模型，多层线性模型等等.对于时间序列来说，需要的时间数据多些效果才会比较好，一般来说至少需要50个，而我们这里只提供6个时间点的数据，显然不合适.后面几种方法都是基于回归的思想发展起来的方法，本文根据数据的特点，将联合多项分布与logistic回归结合起来，来预测产品在每个时间区间的失效率，进而来估计产品的不良返回数.实例的数据结构可以简化为下表的形式：一类月二类月三类月四类月五类月六类月生产数 一月线返回数 三月线返回数六月线返回数十二月线返回数二十四月线返回数三十六月线返回数 表1表1中一类月的数据类型就和07年的各月的数据类型一样的，知道各个月线的不良返回数情况，二类月与08年各月的数据类型一样的，仅仅不知道三十六月线的返回情况，同理知，三类月与09月的一样，四类月与10年一月至六月一样，五类月与10年七、八、九月一样，六类月与10年十月和十一月一样.表中生产数表示这个月生产产品的数目，各个月线的返回数与数据来源中的介绍是一样的.表中表示第i类月生产的产品数，表示第i类月生产的产品数在j月线的不良产品返回数，而作为厂家来说，就希望将这个表中那些空白的部分预测出来.四、模型一：联合多项分布概率模型 1、建立模型和模型求解假设所有产品的生产工艺和技术都是一样的，对于一个产品是否失效可以看作一个随机事件，根据题中数据的分布情况，我们可以假设每件产品在出厂后三个月内月损坏的概率为，第四个月到和五个月之间损坏的概率为，第六个月至第八个月损坏的概率为，第九个月至第十四个月之间损坏的概率为，第十五个月至第二十六个月之间损坏的概率为，第二十七个月至第三十八个月之间损坏的概率为.那么对于一类月的数据的分布情况，每件产品在一月线时间内损坏的概率为，第四个月到和五个月之间损坏的概率为（一月线内没有损坏，三月线内损坏的概率为），第六个月至第八个月损坏的概率为（三月线内没有损坏，六月线内损坏的概率为），第九个月至第十四个月之间损坏的概率为（六月线内没有损坏，十二月线内损坏的概率为），第十五个月至第二十六个月之间损坏的概率为（十二月线内没有损坏，二十四月线内损坏的概率为），第二十七个月至第三十八个月之间损坏的概率为（二十四月线内没有损坏，三十六月线内损坏的概率为），三十六月线以内不坏的概率为.同理我们也可知道二类月、三类月、四类月、五类月、六类月的概率分布情况.我们可以把表1换成一类月二类月三类月四类月五类月六类月生产数 一月线返回数 三月线增加的返回数六月线增加的返回数十二月线增加的返回数二十四月线增加的返回数三十六月线增加的返回数 表2其中，这样做之后相当于将原始数据中的累积不良返回数换成了各个时间段增加的不良返回数.这样的话，我们可以把一类月的生产数与各月线增加的返回的情况看做一个总数为七项分布，并且一类月的数据的似然函数其中 同理我们知道二类月的数据可以看做一个六项分布，相应的似然函数为而 类似地一直到六类月，分别是五项分布，四项分布，三项分布，二项分布，我们也可以得到相应的似然函数.利用上表中所有数据建立的似然函数为这里面有六个参数，对上述似然函数取对数，然后分别对求偏导可以得到的的极大似然估计的递推关系式为单独用第i类月的数据对的似然估计.当很大的时候，并且很小的时候，的估计值可以近似简化为 此时的估计可以看做是各个月估计的加权平均值，相比较前面递推关系式，可以简便很多计算，而且估计也不会差很远.并且由于一般都很大，比较小，由多项分布的性质知，几乎趋于0,也就是说与可以看做是不相关的，也即是相互独立的.一旦把估计出来，我们就可估计表中空着的数据， 比如二类月的三十六月线增加的返回数，因此二月的三十六月线的返回数的估计就为.2、区间估计我们知道点估计对样本的依赖性很大，而且有很大的误差，并且没办法衡量偏差程度.区间估计则可以按一定的可靠性对待估参数给出一个区间范围，因此我们有必要估出不良返回数的置信区间.多项分布的边际分布是二项分布，所以当和的个数比较多时，我们可以近似认为,同理我们也知道，因此，所以对于概率所以，因此的99%置信区间为 查表知.，但是为了方便我们就取来得到置信度更高的区间.例如的置信度为99%的置信区间为： 因此二类月中三十六月线不良返回数的99%置信区间为 3、模型结果模型一中我们并没有假定分布情况，对的估计完全依赖于样本，借助于47个月的部分数据情况建立的似然函数取对数，再分别求偏导数，借助于excel里面的菜单操作计算功能可以得到估计值, ,.而后来我们得到了三个月之后的一些不良返回数据，其实就是2010年十月和十一月的三月线，2010年七月至九月的六月线，2010年一月至三月的十二月线，2009年一月至三月的二十四月线，2008年一月至三月的三十六月线，我们刚好用这些数据来验证的我们的模型，得到如下表的结果月份及月线估计值置信区间真实值及偏差 08年一月三十六月线53356（52045，54691）60407，-703908年二月三十六月线49346（48074，50617）55072，-572608年三月三十六月线77703（761181，79226）91833，-1413009年一月二十四月线34909（33826，35992）38119，-321009年二月二十四月线42501（41245,43757）36450, 605109年三月二十四月线52178（50791，53565）44800， 737810年一月十二月线25066（24131，26001）24142， 92410年二月十二月线21012（20148，21876）20172， 84010年三月十二月线25120（24256，25984）24550， 57010年七月六月线6291（5782,6800）5323， 96810年八月六月线10361（9710，11012）9905， 45610年九月六月线9351（8711,9991）8316， 103510年十月三月线3761（3422,4100）2672， 108910年十一月三月线6104（5683，6525）6402， -298 表3注：表中的置信区间是99%置信区间，偏差=估计值真实值.有阴影月份的月线表示真实值落在我们预测的置信区间内. 从表3的验证结果可以看出，我们的预测效果并不怎么好，14个月份的月线预测结果只有5次落在置信区间内，准确率不到40%，并且不在置信区间内的月份的偏差大多数都比较大，结果比较糟糕.但是我们观察原始数据中的累积不良率AFR（）发现，各个月在同一个时间点的AFR差异还是很大的，譬如一月线的AFR,07年2月的一月线AFR才0.12%,而08年5月的一月线AFR却为0.53%，可见AFR在一月线这个时间点的时间点的跨度是很大的，再看三月线的AFR，最小的才0.35%，最大的达到1.1%，后面的月线也出现类似的结果，而我们的模型中在产品的失效率只与时间段或点有关，其实也就是只与月线有关，只要是同一月线我们认为失效率都是一样的.而样本数据告诉我们，同一月线在不同的月份的失效率差异很大，所以我们的模型单用纵向时间来来估计失效率是不好的，因此我们还要寻找其他的因素来估计在某个时间段的失效率，因而我们就提出了模型二.五、模型二：基于联合多项分布的logistic概率模型1、模型建立由模型一知，我们数据中失效率不仅仅与纵向时间有关.虽然从可靠性理论来说，一件产品的失效主要是与使用时间有关，但是也不能排除某些协变量会对我们的失效率有影响.由于在某一个月份时对的估计我们是用对应的不良返回数除以该月的生产总数，所以我们就有必要考虑生产总数对的影响，所以就把生产总数当作协变量进我们的模型中.由于的范围都在01，而生产总数的值与他们之间的相差很大，因此我们考虑logistic模型，对作logistic变换，从而可以使得它们的范围变成负无穷到正无穷.对于我们要研究的数据，有K=47个月的生产的产品要来估计不良返回数.假定每个月的观察的时间.由给定的六个月线的时间关系知，各个月累积的不良返回数和新增的不良数可以表示为，其中表示07年一月，表示07年二月，以次类推，表示10年十一月，表示第k月在时刻累积的不良返回数，表示第k月生产的产品在时间段期间即将失效的产品个数，用表示失效率，可以解释为第k月生产的一件产品在时间段的失效率，同时我们用表示第k月生产的产品总数，那么第k月在时刻还没有失效的产品总数可以表示为，根据模型一中我们对数据的纵向分析易知，在时间段期间第k月中将要失效的产品服从二项分布，即. 我们建立logit链接函数 (1)这里是具有未知参数的一个已知函数，是模型的协变量，可能的向量取值，是对于于协变量的回归系数向量.而是时间区间内第k月产品失效率的随机效应.设其中 假定这些随机效应在同一个月中相关，而在任两个不同的月中相互独立.设，这里是连续的时间长度的平方根.我们用m维多元正态分布来建立每个月的相关结构，满足 这里是矩阵，定义为是未知参数确定的时间函数，的选定依赖于这类失效数据的反应时间模型.最简单的形式是线性，即，没有常数项的是为了保证模型的可识别性.对的选择还有很多其它的形式.如二次式，甚至更复杂的.但是并不见得越复杂越好，比如Waterloo大学生物系进行的一系列的大种系动物的毒理实验中多元数据分析来看，他们用贝叶斯的方法来选择时发现，最佳的模型是. 对于此时模型的区间估计，我们的思路和模型一的原理是一样的.这里也从多项分布和二项分布之间的关系出发，结合正态分布来寻找的置信区间.由于不同月份生产的产品在同一个时间段内的失效率是不同的，所以模型一区间估计中的要换成.2、模型结果我们同样利用模型一的数据来拟合模型，并且也同样用那14个返回数据来验证模型，由于此时的数据不多，我们只加入各月的生产总数作为协变量，我们就选择了关于时间的一次函数，利用sas编程得到的结果如下：从方差分析的p值小于0.0001知，模型拟合数据时显著性有效的，而对于未知参数估计的p值也都显著性不为0，所以我们的模型为 (2)也即是 (3)这里表示第k月产品的在时间段的失效率，表示第k月的生产的产品总数.而代表的月线，也就是分布为一月线、三月线、六月线、十二月线、二十四月线、三十六月线，但是在我们的模型中的取值分别为1,、2、3、4、5、6.我们也将分别取3、5、8、14、26、38时的取值来拟合模型，但是最后的预测效果没有我们的模型好.为了直观地评估模型的拟合度，模型的学生化残差的散点图如下： 图1从图1我们可以看出大多数标准化残差点都在两倍标准差之内，只有少数点落在之外，此时对于十四个验证的数据的估计如下：月份及月线估计值置信区间真实值及偏差 08年一月三十六月线62323（60868，63778）60407，191608年二月三十六月线57096（55702，58489）55072，202408年三月三十六月线92920（91151，94690）91833，108709年一月二十四月线31554（30539，32568）38119，-656509年二月二十四月线36815（35682,37948）36450, 36509年三月二十四月线47002（45723，48282）44800， 220210年一月十二月线23437（22528，24346）24142， -70510年二月十二月线19143（18321，19965）20172， -102910年三月十二月线23529（22618，24440）24550，-102110年七月六月线5689（5238,6139）5323， 36610年八月六月线10051（9453，10650）9905， 14610年九月六月线8992（8426,9558）8316， 67610年十月三月线3355（3009,3702）2672， 68310年十一月三月线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年一月的数据变异性较大.但总体来说，比模型一的结果还是好了很多.六、模型改进和解释 1、模型改进由于后来我们得到了这个产品（MB）一些月份三十六个月线的返回数据，所以此时我们就用了每个月一月线、二月线、一直到三十六月的数据来拟合模型二.由于这时的数据多了很多，所以此时的协变量只取每个月的生产总数显然不合适.我们就想是否具有季节性，因此在协变量中加入了月份，但是月份的影响是周期性的，故构造了周期为12的三角函数对失效率的影响，模型为： （4）此时的其中，另外我们还加入了生产总数与月线的交互项，拟合数据得到的结果如下：从sas运行的结果可以看出，模型的检验和系数的检验都可以通过.其中sn表示的是，cn表示的是，故此时的模型可写为：也即是 这里表示第k月产品的在时间段的失效率，表示第k月的生产的产品总数.而代表的月线，也就是分别为一月线、二月线、三月线、三十六月线，但是在我们的模型中的取值分别为1,、2、3、36.而s表示的是将k对12取余，即，得到的标准化残差图： 图2从标准化的残差图2可以看出，数据几乎都落在2倍标准差范围之内，只有极少落在之外.我们后来也做了将标准差之外点剔除再来拟合模型，但是由于偏差数据太少的缘故，对我们结果的影响很小.月份及月线估计值置信区间真实值及偏差 08年一月三十六月线60720（60541，60900）60407，31308年二月三十六月线55229（55069，55389）55072，15708年三月三十六月线92128（91893，92362）91833，29509年一月二十四月线38171（37993，38350）38119，5209年二月二十四月线36873（36646,37100）36450, 42309年三月二十四月线45436（45169，45702）44800， 63610年一月十二月线24738（24469，25006）24142， 69610年二月十二月线20189（19957，20420）20172， 1710年三月十二月线24637（24416，24858）24550，8710年七月六月线5411（5250,5572）5323， 8810年八月六月线9489（9264，9714）9905， -41610年九月六月线8316（8297,8759）8316， 21210年十月三月线3156（2944,3369）2672， 48410年十一月三月线6446（6170，6722）6402，44 表5从表5可以看出，此时由于我们月份之间的跨度比较小，所以此时置信区间的长度也比较短，所以此时的偏差总体上也比之前小了很多.08年一月、三月的三十六月线和09年二月二十四月线的真实值与置信区间的界限最近也就相差一百多，经过小小调整就可以了.从上面的偏差知，偏差最大的也不过六百多，离落入置信区间最大也不过三百左右，如果对精度不是要求很高的话，完全可以通过增大区间的长度来落入置信区间内.2、模型解释 对于上面式（5）中关于的关系式，当变大时，由于前面的系数是负的，也将会变大.这很容易理解，因为当一个月的生产总数增加了，而此时工人的生产任务就加大了，而工人的精力有限，从而生产效率降低，因而就会有更多的次品.当变大时，由于前面的系数是正的，反而会减小.从可靠性理论的角度来说，一个产品的失效的概率分布通常可以看做一个指数函数，由于指数函数的分布函数，我们可假设一个相隔一个时间间隔的概率函数，而求导得所以是一个关于t的减函数，也就是说随着时间的增大，相等的时间间隔内失效的概率是减小的，这刚好也我们的模型是吻合的.而对于月份的影响也是很容易解释的，由于每个月之间的气候、湿度等等自然因素的不同，对产品的保存也会造成一定的影响，故月份也是有影响的，只是说相对来说小点而已.七、两个模型的讨论与分析两个模型从预测的方式来说可以说方式是一样的，都是先预测一段时间内的失效率，再用失效率乘以生产总数得到增加的不良返回数，但是预测失效率的方法是不一样的，并且引入的参数也不一样.模型一中我们的模型假设是在每一件生产产品在出厂后的相同的一段时间内的概率是相同的，只依赖产品的出厂时间，从而借助多项分布来构造样本的似然函数来估计各个区间的失效率.但是例子验证的结果并不理想，区间准确率约为36%，而其他落在区间外的月份的月线的预测都离区间比较远，所以效果并不好.而模型二就引入了协变量来估计失效率.除了考虑出厂时间，我们还考虑了生产总数的影响和月份的季节性影响，甚至还考虑他们的之间的交互，准确率大约维持在50%多，但是我们发现还有几个月的真实值离我们的预测区间只有100左右的偏差，我们的预测区间本来是就比较短，只要区间稍微放大，我们就可以使得准确率超过80%.对于改进的模型来说，偏差最大也就600左右，也就是说如果我们把置信区间的长度放大到1200多，那我们就有可能全部预测对.此时的预测都是基于上一月线已知的不良数加上新增的不良数来预测下一月线的不良数，因此我们就想能不能跳跃性地预测，预测之后几个月线的.因此我们用改进的模型基于二十四月线的数据来直接预测三十六月线的数据，也就是我们要求二十四月线至三十六月线之间增加的不良数，我们利用模型分别算出了二十四月线至二十五月线的增加数，二十五月线至二十六月线的增加数，三十五月线至三十六月线的增加数，再加这十二个间隔的增加数加起来预测二十四至三十六月线的增加数，我们计算出来的结果发现全部偏差很大，甚至还比不上我们之前就只用一、三、六、十二、二十四、三十六月线的数据拟合的模型，也就是说我们这个模型只适合预测一个时间点的数据，并且这个数据只与上一个时间点数据有关的数据类型，与更前面的的数据关联不大，因此也不适合作跳跃性的时间预测.对于模型二中存在小部分数据的学生化残差不在两倍标准差以内的数据，我们剔除了这些点再来拟合模型，结果发现对最后的结果影响不大，也就那些数据对模型的贡献不怎么大.对于模型中出现的生产总数与月线的交叉项，我们也用过不含交叉项的模型去拟合数据，结果和我们上面的差不多，都只是有小小的变动，所以也可以考虑用不含交叉项的模型拟合.从实例的结果来看，尽管模型二的结果相对于模型一的结果来说有所改进，但是结果还不是特别满意，也只能说还可以，虽然改进的模型中预测的最大偏差不过600多，但是我们仔细观察数据，发现从上一月线到下一月线增加的不良返回数也不过就是2000左右，这样来说600也还是比较大的偏差.对于此例的偏差，我们分析可能有以下两方面的原因：（1）一方面可能来自于我们的模型，可能是我们的模型不是太适合这个数据，可能是有某些重要的协变量我们没有找到，这可能涉及到电子产品制造的某些专业知识，也正是我们所欠缺的.（2）另一方面可能来自于数据，我们的数据来源于回收的不良产品数，而我们的模型是应用于产品的失效前提下，这两者是不等同的.因为有些产品失效了，但是可能却没有及时回收.有些人怕麻烦或者比较忙没有及时将自己的失效的产品返厂，可能就不返厂了或者拖到下个月、下下个月甚至更晚返厂.这样的话就会导致数据会有偏差.八、模型的不足与改进本文结合联合多项分布和logistic回归提供了一种解决多个联合不完全纵向数据的新的方法，虽然从实例来说不是特别理想，但是也不失为这种数据类型可以选择的一种分析方法.我们是先通过预测失效率，再来预测产品失效的个数，由于失效率比较小，相应的标准差也就小，所以我们求得的置信区间的长度不是很大，也就说更强调精度，应用价值更好.我们两个模型的区间估计几乎是用一样的方法，模型一中的模型假设某一段时间的实效率是与生产总数是无关的，而模型二中失效率是与生产总数是有关系的，而考虑失效率一般都很小，即使有生产总数的影响，偏差也不会大得很多，本来我们的区间就不长，所以我们认为是可以接受的，当然从假设模型的分布得到模型二的区间估计是更有说服性. 我们的模型二中是借助于logistic模型来估计失效率，也许我们还可以借助于多项式曲线回归、Gompertz曲线回归、二项型指数曲线回归、三项型指数曲线回归等等更复杂的模型来拟合失效率，或许可以得到一些更好的模型.另外，从改进的模型的解释部分可知，我们的模型一定程度可以用指数分布来解释，然而我们知道在产品寿命研究还有一个更复杂点的威布尔分布，如果我们能构造一个模型最终可以用威布尔分布来解释，可能会有意想不到的结果.对于协变量的选取，我们模型选了生产总数和月份，可能还有某些重要的协变量我们没有找到，结合电子制造和工程学的相关知识对我们寻找协变量可能也会有很多的帮助，由于得到数据的时间相对来说比较短，我们没有时间去学习这方面的知识，但应该是一个可以考虑的方面.参考文献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.附件程序一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,i)/a(1,i);end; do i=3 to 7;do j=1 to 47;if a(i,j)=0 then b(i,j)=(a(i,j)-a(i-1,j)/a(1,j);end;end;do i=2 to 7;do j=1 to 47;if b(i,j)=. then do;n=b(1,j);t=i-1;p=b(i,j);y=log(p/(1-p);end;else do;n=.;t=.;p=.;y=.;end;output;end;end;keep n t p y;cards;631301601312841863781089681294977027116827210582811008180120693873512576917610982301013610145386211978261457726162816418779561449110140161410884541129930812878884558118954714507631201461971568121747415720711758762143069297355493419812701191139780973097973097912274802259652224592574969157937056686385991174146672426881312752205017361728205929521505169928423403596943117679725554434011492544872867318235043778369133174011371557505878435226833332296029132135262812931403157111381764164311022264252721225972369425856639555550885634746940444343689374191422810486159851542212380785187748554637765306777776394549582892675791282312987101245266646586507372596269474711529345522989496041290066524618132107578598512413116601023910628135467626844113971147112625818238261772386420484136071497714673109631095411989144141744115679133601198923806233811661010157115601556914264117891414597311063582290000015337949920839143191038120982214531899519848271141536615154249982548041547292904049841379388502869331084290192185120616202742282028176244392069019710409594019229292180782129026739000000000002824516016294922314815120310203480031400311224956625391259514443641102658795203867598761296936461151645165534935867389550000000000000000000000036823202643458128525178053640041521378033794466206331093586400000000000000000000000000000000000;proc reg data=chengxu1;var y n t;model y=n t/p r ;print cli clm;plot student.*Predicted.;run;程序二data chengxu2;array a37,50;array b37,50;pie=3.1415926;do i=1 to 37;do j=1 to 50;input x;a(i,j)=x;end;end;do i=1 to 50;b(1,i)=a(1,i);b(2,i)=a(2,i)/a(1,i);end; do i=3 to 37;do j=1 to 50;if a(i,j)=0 then b(i,j)=(a(i,j)-a(i-1,j)/a(1,j);end;end;do i=2 to 37;do j=1 to 50;if b(i,j)=. then do;n=b(1,j);t=i-1;p=b(i,j);y=log(p/(1-p);m=mod(j,12);nt=n*t;sn=sin(pie/6)*m);cn=cos(pie/6)*m);end;else do;n=.;t=.;p=.;y=.;m=.;nt=.;sn=.;cn=.;end;output;end;end;keep n t p y m nt sn cn;cards;631301 601312 841863 781089 681294 977027 1168272 1058281 1008180 1206938 735125 769176 1098230 1013610 1453862 1197826 1457726 1628164 1877956 1449110 1401614 1088454 1129930 812878 884558 1189547 1450763 1201461 971568 1217474 1572071 1758762 1430692 973554 934198 1270119 1139780 973097 1157973 912274 802259 652224 592574 969157 937056 686385 991174 1094923 970466 737902 1466 724 2688 1312 752 2050 1736 1728 2059 2952 1505 1682 2124 2660 6058 4366 7765 7245 5484 4074 4993 4469 2852 3182 3505 3781 3693 3319 4013 3716 5751 5887 4352 2685 3332 2960 2913 2135 2629 1293 1403 1571 1143 1782 1650 1123 2381 2475 2521 1756 1884 1387 4127 2412 1749 3778 3604 3280 4082 5436 2679 2753 4204 4426 9482 7675 11683 11013 9487 6019 7454 6150 4665 5024 5295 5703 6427 6129 6809 5933 9730 9219 7573 4254 4203 6285 4851 4067 4615 2553 3203 2923 2019 3452 3073 1897 4321 4652 4535 0 2527 2122 5972 3694 2585 6639 5555 5088 5634 7464 3411 3867 6109 6731 14447 10631 16234 15426 12557 8038 8918 8538 6363 6530 6784 7768 9459 9584 8932 7588 12828 12996 10124 5268 6465 8650 7372 5962 6948 4711 5296 4564 3000 4982 4173 2672 6402 8154 0 0 3723 2908 7728 4940 3731 8512 7515 6804 7466 9367 4453 4780 8174 9236 19178 13599 20484 18691 15740 9565 10900 10659 8013 7592 8207 10141 12487 12091 10553 9181 17785 17050 11783 7352 8237 11060 9737 7804 9487 6372 7257 6055 4026 6377 5400 3796 9038 0 0 0 5373 3753 11705 6244 4765 10685 9447 8558 8981 11441 5692 6169 11231 11640 23244 16128 23776 21614 17894 11640 13061 12792 9223 9054 10022 12329 15382 13962 11971 10739 22022 19150 14813 8752 9932 13227 11781 9904 11995 8112 9190 7362 4687 7877 7047 4754 0 0 0 0 6652 4618 13210 7578 5986 12413 11660 10221 10556 13394 6889 7971 13299 14176 26650 18510 26620 23937 20782 13854 15256 14672 10951 10954 12000 14422 17452 15685 13368 11998 23811 23390 16610 10159 11560 15570 14264 11796 14170 9766 10681 8301 5323 9905 8