分位数回归在大气环境中的应用

1、分位数回归在大气环境中的应用摘要：文中给出了分位数回归的方法，并将其应用于大气环境浓度中进行了实证分析。文中采用纽约市1973年5月到7月三个月的数据对太阳辐射、温度以及风速对臭氧浓度的影响进行了相关分析，利用用分位数回归模型进行了求解，得到了多种模型结果。关键字：臭氧，臭氧浓度，辐射，分位数回归现今，大气中臭氧层对地球生物的保护作用已经广为人知，它能够吸收太阳释放出来的绝大部分紫外线，使得动植物免遭其危害。但是由于工业的发展而带来的环境污染等因素，使得地球外围的臭氧层日益稀薄甚至出现空洞。为了应对臭氧层逐渐遭受破坏，并由此产生的对动植物的广泛而严重的威胁局面，国际出台多种政策，推广专门针对臭

2、氧层保护的技术和标准，以减少进一步破坏甚至弥补日渐稀薄的臭氧层乃至臭氧层空洞。1987年通过，并于1989年生效的蒙特利尔议定书标志着世界各国共同应对臭氧层破坏的开始，此后，世界各国积极采用先进技术，比如推广使用无氟制冷剂等措施，以减少氟利昂等物质对臭氧层的破坏。世界上还为此专门设立国际保护臭氧层日。然而，大量的保护臭氧层保护环境的宣传给人的印象似乎是受到保护的臭氧应该越多越好，其实并不是这样，如果大气中的臭氧，尤其是地面附近的大气中的臭氧聚集过多，臭氧浓度过高，对人类来说反而是个祸害。臭氧是地球大气中一种微量气体，它是由于大气中氧分子受太阳辐射分解成氧原子后，氧原子又与周围的氧分子结合而形成

3、的，含有3个氧原子。大气中90以上的臭氧存在于大气层的上部或平流层，离地面有1050千米，这才是需要人类保护的大气臭氧层。还有少部分的臭氧分子徘徊在近地面，仍能对阻挡紫外线有一定作用。但是，近年发现地面附近大气中的臭氧浓度有快速增高的趋势，也引起了环境学家的注意。地面附近越来越多的臭氧，同铅污染、硫化物等一样，来源于人类活动，汽车、燃料、石化等是臭氧的重要污染源。在车水马龙的街上行走，我们常常看到有些空气略带浅棕色，又有一股辛辣刺激的气味，这就是我们经常提到的光化学烟雾。臭氧就是光化学烟雾的主要成分，它不是直接被排放的，而是转化而成的，比如汽车排放的氮氧化物，只要在阳光辐射及适合的气象条件下就

4、可以生成臭氧。随着汽车和工业排放的增加，地面臭氧污染在欧洲、北美、日本以及我国的许多城市中成为普遍现象。根据专家目前所掌握的资料估计，到2005年，近地面大气臭氧层将成为影响我国华北地区空气质量的主要污染物。研究表明，空气中臭氧浓度在0.012ppm（百万分之一）水平时这也是许多城市中典型的水平，能导致人皮肤刺痒，眼睛、鼻咽、呼吸道受刺激，肺功能受影响，引起咳嗽、气短和胸痛等症状；空气中臭氧水平提高到0.05ppm，入院就医人数平均上升710。原因就在于，作为强氧化剂，臭氧几乎能与任何生物组织反应。当臭氧被吸入呼吸道时，就会与呼吸道中的细胞、流体和组织很快反应，导致肺功能减弱和组织损伤。对那些

5、患有气喘病、肺气肿和慢性支气管炎的人来说，臭氧的危害更为明显。本文将通过分析太阳辐射、风速和日最高温度与臭氧浓度的关系，用统计学方法从数据中挖掘出可靠的信息，以便在生活中预测户外的地面臭氧浓度，具体而言，就是在一定的温度、日照、风速等气象条件下为人们提供准确的预警信息。一、模型建立1.1分位数回归理论统计学上“回归”一词最初是由高尔顿提出的，高尔顿在研究人类身高遗传的时候发现这么一个规律：父母身高高，其子女的身高有低于父母身高的趋势，父母身高矮，其子女的身高有高于父母身高的优势。所有身高高、低的人的子女都有“回归”到平均身高的趋势。这样就解释了人类身高没有出现两极分化的现象。1886年，高尔顿

6、正式提出了“回归”的概念，并通过他的学生皮尔逊多年的研究和发展，使得这个生物学领域的概念被推广成为统计方法中的重要概念。传统的线性回归模型描述了因变量的条件分布受到自变量 x 的影响过程。普通最小二乘法是回归分析中经典的、应用最为广为的回归方法，其解释与人们的直观想象相一致，且易于计算，它描述了自变量 x 对于因变量 y 的均值影响。如果模型中的随机扰动项来自均值为零而且同方差的分布，那么回归系数的最小二乘估计为最佳线性无偏估计(blue) ；如果近一步随机扰动项服从正态分布，那么回归系数的最小二乘法或极大似然估计为最小方差无偏估计(mvue)。但是在实际的经济生活中，能够满足上述条件的情况并

7、不多见；例如数据出现尖峰或厚尾的分布、存在显著的异方差等情况，这时的最小二乘法估计将不再具有上述优良性且稳健性非常差。最小二乘回归假定自变量 x 只能影响因变量的条件分布的位置，但不能影响其分布的刻度或形状的任何其他方面。由于最小二乘的估计结果稳健性较差，容易受到异常点或缺失点的影响。即使所遇到的问题满足上述情况，大量的宝贵数据通过最小二乘回归分析方法也只能估计出一条回归曲线，得到的信息还比较有限，很不全面。因此，在使用最小二乘回归的同时，人们也在不断的寻找更好的、能够提供更全面数据信息的回归分析方法。1978年，koenker和bassett首次提出了“分位数回归”的概念。1818年lapl

8、ace曾提出中位数回归理论，他们在此基础上把中位数回归推广到成为一般的分位数回归。与最小二乘回归相比，分位数回归的应用条件更为宽松，挖掘出的数据信息更加丰富，它可以度量回归变量在分布中心的影响，还可以捕捉整个条件分布的特征，度量回归变量在分布上尾和下尾的影响。尤其在假设误差为非正态分布的时候，分位数回归估计量的有效性要强于最小二乘估计量，因此，分位数回归弥补了最小二乘回归的不足。在近二、三十年中，koenker，bassett，powell，chernozhukov等人为分位数理论的发展做出了巨大的贡献。koenker和bassett在1978年把响应变量作为其他变量的线性函数，推到出分位数回

9、归系数的渐进分布，发展了线性分位数回归理论；1982年，他们又研究了分位数回归的异方差稳健性检验和线性假设检验，为分位数回归提供了可靠的保证；1986年，bassett等人又研究了分位数回归强相合性等性质；也在这一年，powell基于删失模型提出了非线性分位数回归；kim和white研究了非线性分位数回归估计量的一致性等性质；1987年，koenker 等人提出了针对线性模型的 l 估计法；同年，koenker等人提出了关于分位数回归的有效算法；buchinsky在 1995 年和 1998 年分别讨论了分位数回归模型渐近协方差矩阵的估计方法，和分位数回归一些最新发展，并应用它分析了美国女性薪

10、水结构的变化情况；2000 年，koenker 和 zhijie xiao解决了分位数回归过程中存在的特定推断问题；同年，kim 和 muller二人研究了两步分位数回归的渐近特性；2001 年，tasche对分位数回归的无偏性进行了研究；2002 年，koenker 等人对线性异方差模型的 l 估计法进行了讨论。此外，chernozhukov 和 han hong提出了关于删失分位数回归的三步评估的研究方法；吴建南和 bretschneider 等用蒙特卡罗(monte carlo)方法随机产生出来的 100 个数据集合用以来比较分位数回归和显著加权分析方法的优劣；kottas 和 krnj

11、ajic又提出分位数回归中的贝叶斯非参数模型；chernozhukov提出了极端分位数的概念，推导出极端分位数的渐近分布及其在经济领域的应用。至此，人们基本建立了分位数回归的理论体系。近 10 多年来，分位数回归在国外得到了迅猛的发展及应用，其研究领域包括经济、 医学、环境科学、生存分析以及动植物学等方面。为了说明分位数回归的有用性，我们特介绍两个分位数回归实证分析的例子。koenker 和 machado 分析了 19651975以及 19751985这两段时间内世界主要国家的经济增长情况，模型选取了13个影响经济增长的自变量，通过分位数回归得出结论:对于起初的单位资本产出这一自变量来说，它

12、的全部回归分位系数基本保持不变，这就意味着对于经济发展迅速与缓慢的国家而言，起初的单位资本产出对于经济增长的影响基本相同；但是教育支出占 gdp的比重以及公共消费占 gdp 的比重这两个自变量对于经济发展缓慢的国家影响更加的强烈。chen使用分位数回归方法深入研究了美国8250名男性的bmi (身体质量指数，一种广泛用于测量偏胖还是偏瘦的指标，bmi=体重/身高2)情况，并得出结论：在 220 岁这一快速成长期中，bmi 非常迅速的增加；在中年期间其值保持比较稳定；60岁以后，bmi的值开始减少。这对于如何保持一健康的身体提供了一种非常有效的措施 ,可以在各个阶段中分别采取相应的控制体重的方法

13、。在分位数回归提出及发展的早期阶段，由于本身计算的复杂性，分位数回归并没能像经典的回归分析那样迅速普及，但是随着分位数回归理论的完善，关于分位数回归的算法也在不断地发展中，现在应用最为广泛的有以下三种：(1)单纯形法，即任选一个顶点，然后沿着可行解围成的多边形的边界搜索，直到找到最优点，这种算法的特点决定其适合样本量不大和变量不多的情形；(2)预处理后内点法，即对回归系数预处理后，再用内点法，它适合样本量大于 10 的情形；(3)内点法，即从可行解围成的多边形一个内点出发，但不出边界，直到找到最优点，它对于处理大样本问题时效率比较高。而由于计算机的辅助功效，分位数回归在近几年逐步成为了经济、医

14、学、教育等领域的常用分析工具。对于分位数理论的应用，国外的研究以及应用的文章相对很多，主要有 barnes 和w.hughes利用分位数回归对跨部门公债市场的回报率进行了分析；buchinsky应用分位数回归研究了 1963 年至 1987 年美国工资结构的变化趋势；bouyé 和salmon应用非线性分位数回归研究了汇率市场中不同汇率之间尾部区域的风险相关性；buhai在回顾分位数回归方法的基础上,研究了它的两个主要应用,分别是持续时期模型和循环结构等式模型；leggett 和 craighead利用分位数回归确定了时间的分布和特定风险驱动的影响；whittaker 等在信用贷款记

15、账方面应用该技术进行分析。然而，国内关于分位数回归的研究较少，且大多数都停留在应用上。如曾昭玲等利用分位数回归模型探讨负债融资及权益融资在不同分位数时,受解释变量影响的差异性,以期对于不同融资程度的公司受资本结构变量及众变量影响的差异性进行深入了解；蔡明璋利用分位数回归模型估计家务工作时间的变异分配；荀鹏程等对中位数回归模型进行了介绍,并将其应用于北京市sars发病预测中；季莘等对百分位数回归进行了介绍，并将其应用于制订正常人群血压参考值的研究。此外，李育安介绍了分位数回归的概念、算法以及在经济领域的应用；吴建南和马伟把分位数回归技术与显著加权法通过估计极端行为参数的能力进行了比较。总之，国内

16、关于分位数回归理论及应用的研究在不断地发展和完善。1.2分位数回归基本概念及定理 2.1分位数回归的定义如果在一组规模为的样本中，第个样本的观测值超过了（）个其他样本观测值，而低于个样本的观测值，那么，我们就说样本的观测值排在全体样本的第分位。样本成为分位点，样本观测值为分位数。对于任意分布函数为的实值随机变量x，任意的满足，有则称是x的分位数，也可记为。当时，称为中位数。我们可以把分位数看做是一个优化问题的解。仿照最小二乘方法，解下式由此得到为样本均值。当给定时，如果把的条件均值记为，即是的线性函数，则可以由估计得到，这就是经典的最小二乘回归问题。同样的，求解可以得到，即为样本的分位数。当是

17、给定的时候，的条件分布函数几座，其反函数为。定义为样本条件分位数函数，其中为维向量。则可以由估计得到。这个思想最早是由koenker和bassett提出的。更深一步，我们可以把上式整理为其中，为系数向量，并且它是的函数，随的变化而变化。2.2 分位数回归的渐进理论对于独立随机变量，其分布函数分别是。而的个观测值组成了向量。假设分位数函数关于协变量为线性的，其中为参数向量。变量的条件分布函数可以写成。定理1 当满足条件（1） 分布函数绝对连续，而且其连续密度函数在点的范围是；（2） 存在正定矩阵，使得 时有，对于假设误差服从独立同分布的模型有，这里。2.3分位数回归的适用情况和优势分位数回归主要

18、针对的对象是多个输入变量、单个输出变量的模型，也就是我们常见的多个自变量，一个因变量的情况。分位数回归同其他大多数的统计分析技术一样，首先需要选定研究模型，确定诸如自变量和因变量相互关系、自变量和自变量相互关系等研究假设。然后，根据原始的测量数据，绘制散点图，观察原始数据点的分布，初步判定数据分布的类型（线性分布或者非线性分布等）。当数据分布中存在较多孤立点或者奇异点的时候，分位数回归的应用意义比较大，效果相对于其他回归分析要较好。在使用分位数回归分析的时候，对于分位数的选取也是很重要的。可以选取特定一个分位数，也可以选定一组或者满足一定变化范围的分位数，例如0.25、0.5、0.75等。分位

19、数选取的个数由样本容量的大小来决定，容量大的样本选取的分位数个数要多一些，容量小的样本可以选取较少个数的分位数。除了变量个数一定要小于样本容量以为，分位数回归分析中并没有对于关于自变量个数上限的具体规定。但是，如果分析当中的变量个数接近样本容量的时候，就会牺牲自由度，函数模型的预测能力会更多的被多余的变量的影响，这就要求当样本容量较小的时候适当的选取自变量。当样本容量较大的时候，就比较不用考虑这个问题。分位数回归是采用最小化加权残差绝对值这和的方法来估计参数的，它的优点主要有以下几点：分位数回归对于模型中出现的随机误差项的分布不做要求，这样就使的分位数回归的适用范围更加广泛；分位数回归具有很好

20、的稳健性，其参数的估计值不受或很少受到异常点的影响；分位数回归的分位点取值不同时，所得到的参数估计值也会不同，也就是取多个值，就可以挖掘数据中相对全面的信息。二、实证分析3.2 数据分析我们选取纽约市1973年5月9月这111天的臭氧浓度（ppb十亿分之一），风速（mph每小时英里数），日最高温度（f华氏温度）和太阳辐射（langley兰利）的日测数据为原始数据。样本容量为444个，对这些数据的基本描述统计量如表3.1所示:表3.1 基本描述统计量臭氧浓度太阳辐射温度风速均值42.0991184.801877.79289.938739方差1107.298308.74290.8203112.66

21、803标准差33.2759791.15239.529973.559218中位数31207799.7峰度系数1.265266-0.4929331-0.22819120.459958偏度系数1.31654-0.9164162-0.6430780.3500279下表是臭氧浓度，太阳辐射，温度和风速原始数据得到的关系系数矩阵，由此可以看出，臭氧浓度与太阳辐射、温度正相关，与风速负相关，而且太阳辐射、温度和风速三个元素之间也并不是相互独立的。表3.2 变量间的相关系数ozoneradiationtemperaturewindzone10.34834170.6985414-0.612951radiatio

22、n0.348341710.2940876-0.127366temperature0.69854140.29408761-0.497146wind-0.612951-0.127366-0.4971461 下面是根据温度和风速的划分下臭氧浓度关于太阳辐射的散点图。图3.1图3.1中o代表臭氧浓度，r代表太阳辐射，w代表风速，t代表温度。本图讲温度和风速划分为四个区间，且为了显著性，每个区间都有重叠。从图中我们可以看出，在温度和风速的划分下，臭氧浓度和太阳辐射的散点图呈现了某种线性趋势。3.3 建立模型经过分析，我们知道，臭氧浓度和太阳辐射呈现某种线性趋势，因此，我们采用下面的现行条件分位数回归模型

23、建模。其中，参数估计值有分位数函数估计得到。表3给出了不同的温度、风速条件下，时，和的估计值。表3.321温度小于80，风速小于102.701298700.051948058.929411760.0470588218.32520325 0.0406504121.55605380.103139062.56451613-0.06989247温度大于80，风速小于102.77777780.1717172-1.00000000.236220513.5291005 0.30687832.53191490.468085118.00000000.4444444温度小于80，风速大于106.276752770

24、.014760158.520.0212.06896552 0.0344827614.116591930.0493273522.513605440.03741497温度大于80，风速大于106.59197320.100334412.810810810.0810810824.92920354 0.0796460245.463768120.0797101444.54545450.0909091从上面的表中可以看出，同一组数据在的取值不同时参数和的估计值是不同的，这比只考虑单一方面性质的线性回归更好更全面的反应了数据各个局部的信息。以温度和风速为条件，将臭氧作为太阳辐射的函数，用分位数回归方法拟合，见

25、图2。图2各图中的五条黑色直线依次是分位数回归直线，红色直线是最小二乘回归直线。图3.2这样我们根据表3，可以得到分位数回归直线方程。例如时，温度小于80f，风速小于10mph，臭氧浓度与太阳辐射的关系是（图2左上图最下面黑色直线）；温度大于80f，风速小于10mph，臭氧浓度与太阳辐射的关系是（图2右上图最下面黑色直线）；温度小于80f，风速大于10mph，臭氧浓度与太阳辐射的关系是（图2左下图最下面黑色直线）；温度大于80f，风速大于10mph，臭氧浓度与太阳辐射的关系是。（图2右下图最下面黑色直线）由此可知，当温度小于80f，风速小于10mph，的时候，太阳辐射每增加100langley

26、，臭氧浓度增加5.19ppb；当温度大于80f，风速小于10mph的时候，如果太阳辐射为200langley的时候，臭氧浓度不超过37.12ppb的概率为10%。同时时，温度小于80f，风速小于10，臭氧浓度与太阳辐射的关系是（图2左上图最上面黑色直线）；温度大于80f，风速小于10mph，臭氧浓度与太阳辐射的关系是（图2右上图最上面黑色直线）；温度小于80f，风速大于10mph，臭氧浓度与太阳辐射的关系是（图2左下图最上面黑色直线）；温度大于80f，风速大于10mph，臭氧浓度与太阳辐射的关系是（图2右下图最上面黑色直线）。这也就是说，当温度大于80f，风速小于10mph，太阳辐射值为300

27、langley时，臭氧浓度不超过151.33ppb的概率为90%，如果臭氧浓度在短期内频繁超过了151.33，也就可以说明当地环境发生了意外变化，需要当地政府和民众引起注意。由上述四个式子以及图2可以看出，当温度较高，风速较小的时候，臭氧浓度正比于太阳辐射上升的速度最快，达到的浓度值相对较高，当温度较低，风速较大的时候，臭氧浓度较低，而且正比于太阳辐射上升的速度也较慢。这也说明，夏季无风的晴天城市地面臭氧浓度会出现较高情况，应该提醒市民出行时避开最高温时段，以减少地面臭氧浓度高而对人们身体健康造成的伤害。而冬天的风天地面臭氧浓度较低，比较不会对户外出行的人们造成身体上的影响。表3.4是我们以温

28、度小于80f，风速小于10mph的原始辐射数据作为检验集，列出不同分为数的分位数回归预测结果、最小二乘回归预测结果与原始臭氧浓度数据的对比。表3.4 两种方法结果的对比辐射原始数据最小二乘臭氧浓度原始数据19012.571417.870526.048741.142549.284932.2911411188.831114.482323.121933.726454.317227.47503629918.233722.999930.479652.394641.666739.58212325616.000020.976428.731747.959644.672036.70581629017.76622

29、2.576430.113851.466342.295738.980111444.987011.000020.113826.094159.489222.52511183.11689.305815.650422.381762.005220.11711254.000010.105819.341424.134560.817221.2542422314.285719.423527.390244.556046.978534.498411527917.194822.058829.666750.318843.064538.224337374.623310.670519.829225.372259.978522

30、.05692073.06499.258818.609722.278062.075220.050216515.350611.329420.398326.816159.000022.9933591158.675314.341123.00033.417054.526827.27432321213.714218.905826.943043.421547.747333.76264525916.155821.117628.853648.269044.462336.90652423815.064920.129428.000046.103145.930135.50182823715.012920.082327

31、.959346.000046.000035.434946143.42859.588218.894323.000061.586020.51842319312.727218.011726.170741.461849.075232.4917301329.506415.094123.650435.067253.408628.344518结语文章用分位数回归对模型进行了求解分析，得到了多种模型结果，相对于单独的最小二乘法更能全面的提取数据中的信息，可以给决策者提供更全面的信息支持。分位数回归对于异常值不敏感，这点也是最小二乘方法所不具备的，因此得出的结果稳健性更高，可信度也较高。参考文献1 koenke

32、r,r.and bassett, g. regression quantiles. econometrical, 1978, 46(1): 33-502 laplace, p. s. theorie analytique des probabilites. editions jacques gabay, paris, 18183 koenker,r.and bassett, g. the asymptotic distribution of the least absolute error estimator. journal of the american statistical assoc

33、iation, 1978, 73: 618-6224 koenker,r.and bassett, g. robust tests for heteroscedasticity based on regression quantiles. econometrica, 1982a, 50: 43-615 koenker,r.and bassett, g. tests of linear hypotheses and l1 estimation. econometrica, 1982b, 50:1577-15846 bassett, g. and koenker, r. strong consis

34、tency of regression quantiles and related empirical processes. econometric theory. 1986, 2: 191-2017 powell, j.l. censored regression quantiles. journal econometrics. 1986, 32: 143-1558 kim tae=hwan and halbert white, estimation, inference and specification testing for possibly misspecified quantile

35、 regression. 2002, ucsd economics dept discussion paper9 koenker, r. and s. portnoy. l-estimation for linear models. journal of the american statistical association. 1987, 82: 851-85710 koenker, r. dorey v. computing regression quantiles. applied statistics, 1987, 36: 383-39311 buchinsky，m. estimati

36、ng the asymptotic govariance matrix for quantile regression models；a monte carlo study. journal of econometrics，1995， 68(2)：303338.12 buchinsky，m. recent advances in quantile regression models. the journal of human resources，1998，1：8812613 buchinsky，m. the dynamics of changes in the female wage dist

37、ribution in the usa：a quantile regression approach.journal of applied econometrics，1998，13：1-30.14 koenker，r. and z. xiao. inference on the quantile regression process. econometrica. 2002, 70: 1583-161215 tae-hwan kim and christophe muller. double-stage quantile regressions. econometrica journal. 20

38、04, 7(1):218-23116 dirk tasche. unbiasedness in least quantile regression,200117 koenker，r. and q. zhao. l-estimation for linear heterscedastic models. journal of nonparametric statistics. 1994, 3:223-23518 chernozhuko, v. and h. hong. 3 step censored quantile regression. journal of the american sta

39、tistical association. 2002, 97:872-89719 jiannan wu，stuart bretschneider and fredrick marc-aurele. estimating models of extreme behavior：a monte carlo comparison between swat and quantile regression.chinese public administration review，2002，1(2)：365387.20 kottas a. and krnjajic m. 2005. bayesian non

40、parametric modeling in quantile regression. technical report 2005-06, ucsc department of applied math and statistics.21 chernozhukov, v. nonparametric extreme regression quantiles.1998,working paper, standford.22 chernozhukov, v. conditional extremes and near extremes: estimation, inference and econ

41、omic applications.2002, ph.d. dissertation, standford23 chernozhukov, v. extremal quantile regression. ann. statist., 2005, 33(2): 806-83924 koenker r, dorey v. a remark on computing regression quantiles. applied statistics，1993(43)：410-414.25 portnoy s, koenker r. the gaussian hare and the laplacia

42、n tortoise：computability of squared-error versus absolute-error estimators. statistical science, 1997 (12):279-300.26 koenker and park. an interior point algorithm for nonlinear quantile regression. journal of econometrics，1996，71:265-283.27 barnes, m. and w. hughes. a quantile regression analysis o

43、f the cross section of stock market returns. working paper, 200228 buchinsky，m. changes in the u.s. wage structure 1963-1987：application of quantile regression. econometrical，1994，62(2):405-458.29 bouyé, e. and m. salmon. dynamic copula quantile regressions and tail area dynamic dependence in f

44、orex markets, 2002, manuscript, financial econometrics research center30 i.sebastian buhai. quantile regression: overview and selected application. unpublished work, 200431 leggett, d. and craighead, s. risk drivers revealed: quantile regression.200032 whittaker, j. and whitehead, c., somers, m. the neglog transformation an quantile regression for the analysis of a large credit scoring database. applied statistics-journal of the royal statistical society series c. 2005, 54:863-87833 曾昭玲,郭乃峰,周时如.