现代统计分析方法及软件应用

现代统计分析方法及软件应用肖天贵引言一统计分析的概念二统计研究的内容趋势变异（化）基本（趋势）平均值（场）（变率）距平（距平场）变化（趋势）增减固有周期（规律）周期律模态随机量不固定周期突变趋势分析、突变检验、周期分析由此可分为三类变率统计特征诊断分析预测方法（统计机动力）各种特征的数据及方法共同交互构成了复杂、庞大的统计方法体系。三本科学习阶段介绍的统计分析与预报方法基本教材：黄嘉佑，气象统计分析与预报方法，气象出版社，2000。四较为流行的统计分析及预测方法1.趋势诊断：滑动平均、累积距平、线性倾向、函数拟合、低通过滤、差分过滤、样条函数法、正交多项式、最优分割法、自适应多分辨滤波……2.突变检验：信噪比SNR、滑动t检验法、Cramer法、Mann-Kendall法、K统计量法、累积距平统计量法、双向回归分析法、Pettitt法、Lepage法……3.周期分析法：功率谱、交叉谱、最大熵谱、奇异谱分析、小波分析、频率域特征分析、Fourier分析、滤波、非线性奇异谱分析……4.场的分析：（针对场变率）场空间结构和时间演变规律：EOF、PCA、REOF、单场CEOF、NLPCA.场波动传播特征：POP、CPOP、CSPOP、HilbertEOFMSSA、WPA、NLMSSA.空间结构及时变规律：SVD、CCA、NLCCA、两个场相互关系RA.波动传播特征：CSVD5.气候变率潜在可预报性分析：方差分析、相关链过程统计检验、气候模式生成的模拟效应检验。6.预测方法一般方法：如回归最优回归、岭回归等灰色预测模型、样条回归……单时间序列：如自回归、均生函数、最优气候均态……多时间序列：相空间投影、相空间相似……非线性预测：门限自回归、时变参数、多层递阶、非线性动态系统、非线性统计的动力模式……相似预测：距离系数、相似系数、逐步相似、秩相关相似预测….神经网络场变率：POP预报模型、线性转置模型LTM、BP-CCA法、最大熵与奇异谱相结合的MEM-SSA法、动力-统计耦合方法……上述方法概括与：•谢明恩、张万诚，云南短期气候预测方法与模型，气象出版社，2000.12•魏凤英，现代气候统计诊断预测技术，气象出版社，1999•吴洪宝、吴蕾，气候变率诊断与预测方法，气象出版社，2005。五、本课程的主要任务：1.简单介绍趋势的突变诊断方法方2.介绍场变率的诊断法以EOF和PCA为基础的诊断方法：如REOF、NLPCA、POP、介SVD、冗余分析等。绍以Fourier为基础的诊断方法：如小波、SSA、NLSSA等。3.简单介绍可预报性及预报方法部分软件的调试（包括简单的编写）软件应用计算分析整个课程学习中的知识准备微积分：多元函数的条件极值基础数学线性代数：m维线性空间、正交基、二次型、矩阵运算、特征值、特征向量等。概率统计：主要统计量、显著性检验、基本理念。数值计算：Fortran、Matlab……软件绘图：Grads、Surfer……几种常用分布：•X2分布•t分布22222212222211222212222....21(),2210,12ninarkkiiiiiiniiXxxxxxEXnVXnYXnYXnXXnXNXnXnNXXnn2()(0,1)()XTtnXNYXnYnF分布2222(,)()()()022(2)()4(2)(4)1(,)arXnFFnmXXnYXmYmmEFmmmnmVFmnmmFmnF第一章趋势诊断与突变检测方法•气候诊断研究就是利用各种手段和工具，对气候及其变化进行诊断，从而得出一个关于气候变化序列的较完整的概念，弄清气候变化的基本规律和物力机制，为短期气候预测提供可靠的信息和物理背景依据。§1.1气候趋势诊断的主要方法一、气候变化及其一般特征气候时间序列：随时间变化的一列气候数据构成了一个气候时间序列，它常表现为离散观测得到随机序列。任何一个气候时间序列Xt都可以看成由以下几个分量所构成：其中，为气候趋势分量；为气候序列的固有周期变化；为循环变化分量；是平稳时间序列分量；是随机扰动项。ttttttXHPCStHtPtCtSt二、滑动平均法•又称移动平均法，是气候趋势拟合诊断技术中最基础的方法，相当于一个低通滤波器，用确定时间序列的平滑值来显示变化趋势，是以一连串部分重迭的序列的平均值组成新序列的一种方法。•对样本量为n的序列x，其滑动平均序列表示为：式中的m为滑动长度。•可以证明，经滑动平均处理后，序列中短于滑动长度的周期大大削弱，显现出变化趋势。111ˆ1,2,,1mjijixxjnmm三、累积距平法累积距平曲线是一种由曲线直观判断气候变化趋势的方法。对于序列x，其某一时刻t的累积距平表示为：11ˆ(1,2,)1ttiiniixxxtnxxn将n个时刻的累积距平值全部算出，即可绘出累积距平曲线进行趋势分析。其中四、线性倾向估计法用表示样本量为n的某一气候变量，用表示所对应的时间，建立与之间的一元线性回归。即用一条直线来表示x与其时间t之间的关系，a为方程的截距，b为斜率，a与b可用最小而乘法计算得出。ixitixixitˆ(1,2,,)iixabtin五、函数拟合法函数拟合法就是建立以时间t及其乘幂项为自变量的多元回归方程来描述趋势分量，故又称为回归分析法。•设气候要素观测序列为，其趋势可用一个时间t的多项式来表示，即根据最小二乘法，由最小可求得系数，这就是以多项式函数拟合趋势。六、低通过滤法把序列的高频分量滤去以便突出长期或气候变化趋势的一种方法。12,,nxxx2012ptpHaatatatpn1()ntttxH012,,,paaaa其过滤后序列为：表示过滤前序列，其权重为，称为滑动区间。权重可以相等，也可以不等，但必须满足所有权重值和等于1。平均的过程实质上是压低小波动（即短周期振动）的影响，而突出长周期波动的作用，因而也起到反映长周期气候因子的综合影响的作用。所取的滑动区间越大，过滤掉的短周期越长，但这种过滤会削弱过滤后振动的振幅和改变其位相。一般用权重分布形状来设计过滤权重。其中，二项系数过滤的权重值的分布遵从二项分布，突出滑动中心的作用。七、差分过滤法将序列x的后期值减去前期值即得差分序列：ˆ()()LiiLxtWxti()xtiiW(21)L()(1)()(2,3,,)dxtxtxttn由于差分是一种高通过滤器，它能过滤出含高频振动的综合序列，那么将原序列减差分序列令得一低频段综合序列。()()()(2,3,,)xLtxtdxttn除上述几种诊断方法外，对气候趋势的诊断方法还有正交多项式法，最优分割法和自适应多分辨滤波法等。八、三次样条函数法三次样条函数拟合实质上是一种分段多项式，它以对给定的时间序列进行分段曲线拟合的方式来反映其本身真实的变化趋势人生管好五张纸，幸福一辈子。•学位一张纸，受益一辈子；•投资一张纸，坚持一辈子；•保险一张纸，安稳一辈子；•房产一张纸，奋斗一辈子；•婚姻一张纸，情暖一辈子。541法则是一种较为流行的收入配置方式——收入50%的资金用于日常的花销，满足现在生活的需要；40%的资金投入到房贷、基金和储蓄中；其余10%的资金则投入到保险账户中。动态管理，贵在坚持。§1.2突变检验的主要方法•引言•（1）突变概念通常意义:变化的过程中的某种不连续现象.可理解为一种质变,一种当量变达到一定限度时发生的质变.气候突变分为简单突变和外力突变.简单突变:天气系统内部不考虑外力影响,有周期变化.如初霜冻,季风爆发等。外力突变:无规律,外力强迫.•(2)研究方法突变理论:以做积分为基础的奇点理论,观察某种系统或过程从一种稳定状态到另一种稳定状态的飞跃。统计理论:统计特性的急剧变化,可考察统计特征值的变化来定义突变.如序列的均值突然变化,方差突然变化,自相关性的突然变化等。目前,均值突变点的检验已有一些方法,但是方差突变和自相关突变尚无较好的检验方法.(3)关于几个统计量•T检验••假设则•又具有可加性UtVn(0,1)UN2()yn2212121212()(,,)xyNnn22121212()(0,1)11xyUNnn22222221(1)()(1)insnsxxn22211221222(1)(1)(2)nsnsvnn•故•一信噪比SNR法(yamamoto法)1986年•含义:以相邻气候阶段的平均值的显著性差异定义为气候突变.•方法:滑动地选择某个参照年,分别计算两个气候阶段内的平均值之差和标准差,定义为信噪比•SNR=•作为相邻两气候阶段差异的统计量，若SNR1.0则认为存在气候突变。若SNR2.0则视为强突变。222212121211221122122121212()()()()(2)1111/(2)2xyxytnnnsnsnsnsnnnnnnnn1212xx•理解:(1)信号与噪音•表示两个气候状态总体平均之差•两个气候状态的差值•表两个气候状态差值的总体标准差,仅反映它们差值的平均振动状况.则•考虑总体的平均值是不变的常数,两边取方差得噪音部分方差•是偏差部分的平均变动程度,则表示信号部分.•表信噪比21xxx21xxx()xxxexxe()xx/()xx22()()xe•(2)显著性检验•原假设气候状态变化且•给出统计量•理解•则据令••故存在关系•若则突变显著.0x12()()xx/(22)txstn221221211121()(,)niiixxxNnnn22211211221212()()112xxtnsnsnnnn121212,,nnIH22222121121122(22)/()()1xxxxxtIHsIHsssxsxIHIH(22)tSNRIHtIHtSNRIH(3)例子若取则相当于,超过的显著性水平表明突变若取则相当于超过的显著性水平，表明基准点出现了强的突变.计算步骤:1.确定基准点前后两段子序列长度,一般取2.连续设置基准点,以滑动方式依次为计算SNR,得到信噪比序列SNRi，i=1,2,….n-2IH-13.若SNRi1.0,认为在i时刻有突变发生SNRi2.0认为在i时刻有强突变发生.1210nnIH1.0SNR3.1622.878tt0.011210nnIH2.0SNR6.324t0.000112nnIH计算结果分析:根据信噪比曲线上的点是否超过1.0或着2.0直线判断序列是否发生过突变或强突变,并确定出发生突变的时间.同时,根据信噪比曲线的变化,分析序列的演变趋势,特别是长期演变趋势.应用案例:魏凤英计算中国北半球,全球的T突变，成都地区T,R计算.•二滑动T检验法(MTT法)•含义:若两段子序列的均值差异超过了一定的显著水平,则认为均值发生了质•变,有突变发生。•方法:对于具有n个样本量的时间序列x，人为设置某一时刻为基准点,基准点•前后两段子序列和的样本为n1和n2,平均值为和方差为•和,定义统计量:•给定显著水平,若则认为发生了突变.1x2x1x2x21s22s121212(2)11wxxttnnsnn221122122wnsnssnnitt•理解:1.t统计量的定义完全与数学书上的一致,容易理解.2.若多个突变点连在一起,则形成了一个突变区域3.子序列时段的选择带有人为性,为避免子序列长度造成的突变点的漂移,可反复变动子序列长度进行比较,提高计算结果的可靠性.•计算步骤:1.确定基准点前后两个子序列的长度,一般相同,2.采用滑