106 系统工程理论与实践 2008年1月 在贝叶斯Logistic回归中假定样本观测值Y=(y1,y2,“,y,)服从参数为p的两点分布，即 ·P(y=1)=p(0),i=1,2,…,n 其中 p9）=ga(@)=,t=1,2,n (5) 并且， 月0+名1月++x13 (8,82,…,8.)r=8= 月。+xu月1+…+x2 人月。+x1月+…+xB 变换函数og(p)=log[p/(1-p)],logt-1为logit的反变换，X为n×k的矩阵，B为k×1的列向量.y:,i =1,2,…,n为取值为零或一的随机样本，x:为任意的实数， 例4失效的密封圈 1986年，“挑战者”号航天飞机在发射升空72秒后爆炸，造成7名宇航员罹难，事故原因怀疑是密封圈 失效，下表为航天飞机在不同温度下发射时密封圈失效的记录. 表1 航天飞机在不同温度下发射时密封圈失效的记录 飞行号 复 9 23 10 13 15 4 3 8 17 失效(=1) 1 1 1 0 0 0 0 0 0 0 0 温度（华氏） .53 57 58 63 66 67 67 67 68 69 70 70 飞行号 11 6 7 16 21 19 22 12 20 18 失效(=1) 1 .1 0 0 0 1 0 0 0 0 0 温度（华氏） 70 70 72 73 75 75 76 76 78 79 81 我们用贝叶斯L0gsc回归模型对表1中的数据做分析.设x表示发射时的温度，p(8)表示在不同温 度时密封圈失效的概率，y=1表示密封圈失效，y=0表示密封圈没有失效， 现在已有一组观测值(x,y:),i=1,2,…,23,按照模型假设 eo 为-B(1,p(8,)》并且p(0)=1+e可 先验分布取为：π(P)=1,似然函数为： L(BI Y) l÷+ (6) 后验分布正比于似然函数 我们先用类似于例1的随机游动采样法进行模拟，取提议函数为N(X,-1,1),构造一条长度为10000 的Markov链.模拟结果见图3. 图中显示估计A。的Markov链几乎没有平稳性可言.同时 8 0.95%的接受比率也使得序列的自相关系数非常高，因而，模拟结 0 果很难令人满意.观察发现回归系数B。与月1高度相关.这启发我 4 200040006000800010000 们如果对回归变量进行变换，可能会得到较好的结果，比如用x- 1.0- x作为新的回归变量设B。,P,为变换后的回归系数，则由：月。+ 0.5 0 1=B。+(x-*),可得： 03、 月。=。-，月= (7) 200040006000800010000 有了（门）式，我们就可以方便地转回到初始的回归系数值上.经过 图3上图为A的轨迹图， 变换后，模拟的结果有了很大的改进.从图象上看，P。,B,的轨迹 下图为B,的轨迹图 万方数据106 系统工程理论与实践 2008年1月 其中 并且 在贝叶斯Logistic回归中假定样本观测值Y=(y。，Y：，…，Y。)服从参数为P的两点分布，即 。P(Y；=1)=P(0；)，i=l，2，…，厅 口． p(0；)=logi￡。1(仇)=南，}=1’2’2…，聘 (0，，0：，…，0。)7=邵= 风+菇ll卢l+…+髫l以 岛+髫2lpl+…+菇2以 风+髫。lpl+…+髫础& (5) 变换函数lo础(p)=log[p／(1一p)]，蛔咖。1为蛔咖的反变换，x为n×．|}的矩阵，卢为J}X 1的列向量．Y；，f =l，2，…，n为取值为零或一的随机样本，茗#为任意的实数． 例4失效的密封圈 1986年，“挑战者”号航天飞机在发射升空72秒后爆炸，造成7名宇航员罹难，事故原因怀疑是密封圈 失效．下表为航天飞机在不同温度下发射时密封圈失效的记录． 表1 航天飞机在不同温度下发射时密封圈失效的记录 飞行号 14 9 23 10 l 5 13 15 4 3 8 17 失效(=1) l l l 1 0 O O O O O O 0 温度(华氏) 53 57 58 63 66 67 67 67 68 69 70 ．70 飞行号 2 11 6 7 16 21 19 22 12 20 18 失效(：1) 1 ．1 0 O 0 l 0 0 O 0 0 温度(华氏) 70 70 72 73 75 75 76 76 78 79 8l 我们用贝叶斯Logistic回归模型对表1中的数据做分析．设菇表示发射时的温度，P(口)表示在不同温 度时密封圈失效的概率，Y=1表示密封圈失效，Y=0表示密封圈没有失效． 现在已有一组观测值(戤，Y；)，i：1，2，…，23，按照模型假设 Y；一B(1，p(B))并且p(∞=亡i≥ 先验分布取为：7t"(卢)毫1，似然函数为： 郴!叭熏(燕斋)^(d‰r ㈤ 后验分布正比于似然函数． 我们先用类似于例1的随机游动采样法进行模拟，取提议函数为Ⅳ(五一。，1)，构造一条长度为10000 的Markov链．模拟结果见图3． 图中显示估计岛的Markov链几乎没有平稳性可言．同时 0．95％的接受比率也使得序列的自相关系数非常高，因而，模拟结 果很难令人满意．观察发现回归系数凤与角高度相关．这启发我 们如果对回归变量进行变换，可能会得到较好的结果，比如用茗一 二作为新的回归变量．设∥。，∥。为变换后的回归系数，则由：风+ 筇l=∥o+(菇一菇)∥l可得： 母o=F q—xF l，pl=F l q、 有了(7)式，我们就可以方便地转回到初始的回归系数值上．经过 变换后，模拟的结果有了很大的改进．从图象上看，风，p。的轨迹 0 2000 4000 6000 8000 10000 图3上图为岛的轨迹图， 下图为风的轨迹图 o 4 o'r m帖o” 万方数据