地理数据统计分析与建模2014

地理数据分析程立cllgg1@QQ:4977442492第一章绪论为什么地理研究中需要数学方法？在地理研究中如何应用数学方法？本课程需要掌握哪些数学方法？3为什么地理研究中需要数学方法？地理学发展历史的必然数学方法在地理学研究中起到重要作用4地理学的历史与数学工具的运用古代地理学——几何学(Geometry测量大地) Geo- Geography 地理学

Geomorphology地貌学

Geobotany 地植物学

Geoinformatics：？

5Geoinformatics–whatisit?

-definedasthescience,technologyandartofcollecting,storing,andanalyzinginformationabouttheearth’ssurface(ocean,land,environment,peopleandnaturalresources)anddisplayingtheanalysedresultsdigitally; -alternativeterm:geographicinformationscience;6地理学的历史与数学工具的运用古代地理学

描写地理事件，地理事实，积累知识(山海经，水经注)

例如：海平面变化的描述7颜真卿

(公元708--784年)

“麻姑自言∶接侍以来，已见东海三为桑田。向到蓬莱，水又浅于往昔会时略半也。岂将复为陵陆乎?方平笑曰∶圣人皆言海中行复扬尘也。”

摘自《麻姑仙坛记》8沈括

(公元1024--1093年)

予奉使河北，遵太行而北。山崖之间，往往衔螺蚌壳及石子如鸟卵者，横亘石壁如带。此乃昔之海滨，今东距海已近千里。所谓大陆者，皆浊泥所淹耳。舜殛鲧于羽山，旧说在东海中，今乃在平陆。

摘自《梦溪笔谈》91832年莱依尔对意大利波佐利大理石柱上的海蚀痕迹的解释10地理学的历史与数学工具的运用近代地理学——统计学 对地理现象进行开始运用数学工具描述现代地理学——多种数学工具 对地理现象进一步定量研究，以揭示地理现象的运动规律与形成机制举例：对海平面描述说明的近现代方法（时间序列分析）112亿5千万年来的全球海面变化过程12最近26万年来的海面变化曲线13最近15万年来的海面变化曲线14

我国棉铃虫与北太平洋海温场月均距平、前期ENSO指标的相关关系，建立棉铃虫的长期预报模型。（线性回归，相关分析）

现代地学研究方法举例

——棉铃虫的长期预警模型15德州郓城丰县南京棉铃虫卵量数据德州,22年(1978~1999)郓城,26年(1974~1999)丰县,20年(1980~1999)16棉铃虫的卵量与前期的海平面温度（SST）距平、ENSO指标存在显著或极显著相关关系不同地区的棉铃虫卵量与前期海温距平存在时空的差异性17山东德州棉铃虫第3代卵量与前期海温相关的时空分布Red:正相关Blue:负相关深色:p<0.01浅色:p<0.0518山东德州棉铃虫第3代卵量与前期海温相关的时空分布19山东郓城

江苏丰县JAN2、FEB2、MAR2均存在大面积的正显著相关区域20

棉铃虫第三代卵量与JAN2北太平洋海温相关的空间分布21Correlationcoefficient棉铃虫卵量与前期各ENSO指标的相关关系22棉铃虫卵量预报模型的建立和验证23OCT2XNOV2X最佳预测模型：Y=-0.0482OCT2–0.6152NOV2历史回检率：70%预测准确率：100%（1994-1999）24Niño-4区的SST距平是预测丰县、郓城和德州三地棉铃虫卵量的最主要因子，其他的ENSO指标的预测能力较差用ENSO指标可以提前15-25个月预测，模型的历史回检率为70%，预测准确率为78%模型预测能力25结论：SST和ENSO指标是害虫大爆发的主要长期预测因子；利用北太平洋海温场SST可以提前20-27个月对棉铃虫作出预测；利用ENSO指标可以提前15-25个月对棉铃虫作出预测。26地理与气象关系的分析同理，我们可以分析其他地理与气象关系重庆高温干旱与三峡工程的关系汶川地震与气象的关系2008年年初的雨雪冰冻天气2010年中国的极端天气2011年西南大旱（旱震关系）日本海啸（3.11）2012年美国大旱2013年云南干旱、雅安地震随机向量简介——多元统计分析的预备知识(现代地理统计)28随机向量的定义若向量ξ=(ξ1,ξ2,...ξn)

中各分量ξi为随机变量,则称该向量为n维随机向量.以随机变量为元素的矩阵称为随机矩阵各个波长的光谱反射率变化构成一个随机向量29随机向量数学期望若

为一n维随机向量,则它的数学期望为

30随机矩阵若矩阵的各个元素为随机变量，则称矩阵X为随机矩阵，31随机矩阵的数学期望

随机矩阵X的数学期望定义为随机矩阵的期望是不是随机矩阵？32随机矩阵数学期望的性质以下AB为普通常数矩阵,X为随机矩阵1．2．3．4．利用定义证明以上公式（作业1）33随机向量的协方差阵(1)若为一n维随机向量,

则:

为该随机向量的协方差阵.34由V(x)的定义可以推出（留作作业2）其中其中D(xi)为随机向量x各分量的(均)方差cov(xixj)为随机变量xi与xj的协方差随机向量的协方差阵(2)35协方差阵的基本性质(1)性质1由于cov(xixj)=cov(xjxi)，所以V(x)是n阶实对称阵性质2cov(xixj)反映了随机变量xi与xj之间的线性相关情况,如果随机变量xi与xj相互独立，则cov(xixj)=0 V(x)成为对角阵

36协方差阵的基本性质(2)性质3：证明：37协方差阵的基本性质(3)

性质4：V(x+c)=V(x)其中c为常数向量（留作作业3）性质5：若y为m维随机向量，x为n维随机向量，V(x)为x的协方差阵,A为m×n阶常数矩阵，并且满足y＝Ax，则：

V(y)＝AV(x)A＇ 证明：

V(y)＝E[(y－E(y))(y－E(y))＇]

＝E[(Ax－E(Ax))(Ax－E(Ax))

＇]

＝E[(A(x－E(x)))(A(x－E(x)))＇]

＝AE(x－E(x))(x－E(x))＇A＇＝AV(x)A＇38随机向量的相关阵定义若为一n维随机向量,则它的相关矩阵定义为:其中：为xi和xj的相关系数

39随机向量相关阵的性质由于rii=1,rij=rji，因此，相关阵为对角线元素为1的实对称阵40随机向量相关阵的性质

若随机向量各分量互不相关则相关阵是一n阶单位阵（rii=1,rij=0）41随机向量相关阵的性质协方差阵与相关阵关系

注意到在相关阵中：其中：V(x)＝DR(x)D

其中：随机向量相关阵展示了一个随机向量内部各随机分量之间的相关性，以及每个随机分量的变动情况（方差）42两个随机向量之间的互协方差阵

定义：若

分别为n维m维随机向量，则向量x和y的互协方差阵为：展开可得到：

43互协方差阵性质对称性：cov(x,y)＝(cov(y,x))＇

Vxy=Vyx＇通常记cov(x,y)为Vxy若随机向量u与x之间满足u=Ax,随机向量v与y之间满足v=By则有：Vuv=AVxyB＇（留作作业4）44互相关阵定义定义：若

分别为n维和m维随机向量则x和y的互相关阵为：

其中r(xi,yj)为随机变量xi与yj

的相关系数,显然:Rxy=Ryx＇多元线性回归——一元线性回归46变量之间的联系变量之间的联系分为确定性关系和非确定性关系(现代统计方法与应用何晓群人大出版社)确定性关系:具有精确的函数关系例如：欧姆定律(V=IR)，自由落体运动公式等非确定性关系:存在密切关系但是没有密切到一个可以完全确定另一个的程度非确定性的原因:存在尚未认识的影响因素试验或测量误差其他种种偶然因素等使得一个或一些随机变量取定值后,其余变量取值带有一定的随机性,不能以确定值与之对应.47变量之间的联系统计关系(相关关系):统计学中把变量间具有密切关联但是不能用函数关系精确表达的关系称为变量间的统计关系或相关关系确定性关系与相关关系之间的转化 由于存在测量误差等原因,确定性关系往往通过相关关系表现出来,许多物理化学定律都是先得到相关关系,在相关关系的启发下逐步发现了这些定律.也就是说当对事物的内部规律了解更加深刻的时候,相关关系可能转化为确定性关系.对于某些统计关系，我们可以利用统计方法建立模型4849坎儿井示意图50回归方程的基本概念回归方法是一种数据挖掘的方法,能从数据中发现有用信息

(是事物的表面规律不是内部机理）数据表格你能从这张表中发现什么规律?年序 最大积雪 灌溉面积y(千亩)

深度x(尺)1 15.2 28.62 10.4 19.33 21.2 40.54 18.6 35.65 26.4 48.96 23.4 45.07 13.5 29.28 16.7 34.19 24.0 46.710 19.1 37.451回归方程的基本概念进行一般加工计算均值(平均最大积雪深度18.5灌溉面积36.53)标准差(最大积雪深度4.773灌溉面积8.746)相关系数只能描述数据的一般特征用于预测还不够好52回归方程的基本概念散点图很容易发现两者存在线性关系,通过计算两者的相关系数也可证明这一点.53回归方程的基本概念通过线性回归可以拟合得到一根最好的趋势线54回归方程应用的另一个例子

直观经验→采集数据→经验验证身高（x）143145146147149150153154155156157158159160162164裤长(y)8885889192939395969897969899100102斜率的意义:整个身高每增高一个单位，腿长贡献其中71.9%5556地理研究的例子5758人口密度与光能辐射596061回归方程的基本概念“回归”（Regression）名称的产生背景来自于英国统计学家F.Galton(1822-1911).Galton在和他的学生K.Pearson(1856-1936)研究父母与子女身高关系的遗传问题时,观察了1078对夫妇,绘制了散点图,计算出一个回归直线方程y=33.73+0.516x

其中,y是成年儿子的身高,x是父母的平均身高.另外：Galton对样本的平均值研究发现： 样本中，父母的平均身高为68英寸（） 儿子的平均身高为69英寸（）似乎表明：人类后代的平均身高总比上一代高1英寸，是这样吗？62进一步研究发现：当双亲的平均身高平均值为72英寸时（大于均值），他们的儿子的平均身高为71英寸（）当双亲的平均身高平均值为64英寸时（小于均值），他们的儿子的平均身高为67英寸（）人类身高出现的所谓的“回归效应”后来人们把Galton所计算得到的直线方程称为回归方程，相应的统计分析称为回归分析实际上，并非所有具有线性关系数据都有回归效应，这样称呼是历史原因造成的回归方程的基本概念63一元回归模型定义：假定（最大积雪深度）x,（灌溉面积）y具有如下线性关系:y=α+βx+εαβ是未知常数,称为回归系数ε表示其他随机因素对y的影响,是一个随机变量.称上式为一元线性回归模型的线性公式变量y称为因变量或响应变量变量x称为自变量或预报变量上式表明：变量y包含有随机成分和线性成分一元线性回归模型是处理两个变量之间关系最简单的模型。虽然简单，但是从中可以了解回归分析方法的基本思想方法和应用。64一元回归模型

实际上，如已经获得n组观察值(x1，y1),(x2，y2),...,(xn，yn)

(比如以上例子中10年的数据） 即自变量x分别取值：x1，x2，...，xn

因变量y分别对应取值：y1，y2，...，yn

符合一元回归模型，则有(xi，yi)满足：

yi=α+βxi+εi65一元回归的前提条件

(Guass-Markov条件)假定1

εi服从正态分布并且E(εi)=0,Var(εi)=σ2(i=1,…,n)现代统计方法与应用何晓群人大出版社假定2

n组数据（样本）是独立观察的,因而,ε1,ε2,ε3,...εn是相互独立的

cov(εi,

εj)＝0（i≠j）注意：模型中假定条件和线性公式一个都不能少66Guass-Markov条件E(εi)=0

表明误差没有任何系统趋势 因而：E(yi)=α+βxi

观测值yi与其期望E(yi)之间的波动εi是完全随机的，没有任何趋势，不是yi的函数，也不是xi的函数Var(εi)=σ2

表明不同次的观测yi在其期望E(yi)附近的波动程度（散布程度）是一样的n组数据（样本）是独立观察的，实际情况中，此条假设较易满足67一元回归模型原理1．已知条件：已经获得n组观察值(x1，y1),(x2，y2),...,(xn，yn),如果它们符合一元线性回归模型，则有：yi=α+βxi+εi并且应有:E(εi)=0,Var(εi)=σ2

通常假定n组数据是独立观察的，因而ε1,ε2,ε3,...εn相互独立的。682．求变量y的数学期望(消除随机误差)对y=α+βx+ε两边求数学期望得到：E(y)=α+βx通常用随机变量y的数学期望E(y)作为y的估计值，记为，故有：即：或同理：对于观测值yi的取值同样有E(yi)=α+βxi或一元回归模型原理69一元回归模型原理3.一元回归的目标：拟合“最佳”直线通过n组观察值来估计α与β，通常用最小二乘法或极大似然法估计通常将，写为：其中记a，b分别为α与β的估计值不同的方法和标准可能得到不同的α与β的估计值。70最小二乘法原理由可知，yi与它的估计值不相等存在一个偏差，称之为残差用ei表示，是观测值与回归值之间的偏差。71最小二乘法原理自然，可用残差平方和用于度量观察值与回归直线之间的接近程度最小二乘法就是使残差平方和Q最小而估计得到的a，b值的方法，就是选择适当的ab使Q达到最小72因为Q(a,b)是关于a，b二次函数，所以Q的最小值是存在的,其中a,b是未知数那么，如何求a，b？最小二乘法原理73最小二乘法原理求极值得到a，b根据微积分中求极值的方法得到方程组：74最小二乘法原理75最小二乘法原理方便起见，记：76最小二乘法原理为求解ab需要证明以下引理：证明：证明：79最小二乘法原理将的第一式化为：即：80最小二乘法原理将代入第二式消去a整理得到：Lxxb=Lxy即b=Lxy/Lxx再由得到a以下是具体的整理过程：代入即b=Lxy/Lxx82最小二乘估计的性质对于一元回归方程：已经获得n组观察值(x1，y1),(x2，y2),...,(xn，yn),如果它们符合一元线性回归模型，则可求得回归方程为：其中b=Lxy/Lxx那么，我们得到的估计值究竟好不好呢？采用何种标准评价？83最小二乘估计的性质无偏估计概念若t是参数T的一个估计，且满足E(t)=T，则称t为T的无偏估计最小二乘估计对α，β的估计具有无偏性：可以证明（多元回归中证明）：a，b是αβ的无偏估计，即：E(a)=α,E(b)=β就是说若用同样的估计方法对α，β作多次估计,a，b的平均值将趋于α，β并且由此可以证明y的估计是E(y)无偏估计84最小二乘估计的性质y估计无偏性的证明表明回归值（估计值）的均值看作实际观察值的平均值注意此性质与Guass-Markov条件的内在联系85最小二乘估计的性质关于a，b的方差（波动性）

由：得到再由yi的独立性及Var(yi)=σ2,得到（注意下标i，j）

86最小二乘估计的性质同样根据可以得到87最小二乘估计的性质由此可以看出：斜率b的方差（随机变量取值波动的大小）不仅与随机误差ε的方差有关而且与x的方差有关截距a的方差与x的方差，随机误差的方差，观察值（样本）个数有关。88作业4作业：利用公式计算例子中的一元回归直线方程(要有计算过程)

最大积雪深度x(尺)灌溉面积y(千亩)115.228.6210.419.3321.240.5418.635.6526.448.9623.445.0713.529.2816.734.1924.046.71019.137.489作业5仿照的推导过程（提示：展开后推导），证明：91回归方程的显著性检验显著性检验的目的任何数据即使没有线性关系都有可能建立回归方程,这样是没有意义的,没有反映变量之间的实际关系.一方面,要建立从经验上认为有意义的方程另一方面,需要用数学方法对方程的显著性进行检验.92数据的相关性93数据的相关性94回归方程的显著性检验

(相关系数法)定义线性相关系数r相关系数与回归系数b的关系(符号相同)

计算例子中的线性相关系数为:r=0.9894,因此,最大积雪深度与灌溉面积有密切的线性关系,并且是正相关.

95回归方程的显著性检验

(相关系数法)相关系数检验表的使用

1．计算得到的相关系数的绝对值必须大于表中的值

2．通常，r的绝对值大于表中α＝0.05的相应值,但是小于α＝0.01的相应值时称x,y有显著的线性关系,

如果r的绝对值大于表中α＝0.01的相应值称有十分显著的线性关系,

小于表中α＝0.05的相应值称x与y没有显著的线性关系.

96

相关系数显著性检验表

P(|γ|>γα)=ααN-20.100.050.020.010.001123456789...0.98770.90000.80540.72930.66940.62150.58220.54940.5214...0.99690.95000.87830.81140.75450.70670.66640.63190.6021...0.99950.98000.93430.88220.83290.78870.74980.71550.6851...0.99990.99000.95870.91720.87450.83430.79770.76460.7348...0.99990.99900.99120.97410.95070.92490.89820.87210.8471...

例子中,n=10,表中α＝0.05的相应值(n-2=8)为0.632,α＝0.01的相应值(n-2=8)为0.765而r=0.9894>0.765,因此,最大积雪深度与灌溉面积有密切的线性关系.97回归方程的显著性检验

(相关系数法)相关系数法的缺点与数据组数n有关组数小时r容易接近1比如:只有两个点的情况（相关系数总为1)98回归方程的显著性检验

(F检验)a.记Lyy总离差平方和为S总

b.回归离差平方和S回c.残差平方和S残99100回归方程的显著性检验(F检验)证明:S总=S回+S残已知：101回归方程的显著性检验(F检验)替换yi的估计根据Q（a，b）中的第一式替换a103回归方程的显著性检验(F检验)公式S总=S回+S残的解释:y的偏差由两个原因造成:一个是x的变化引起y的变化(S回线性成分),一个是由随机误差造成(S残).由此可见S回所占S总的比重越大，回归效果越好104回归方程的显著性检验(F检验)F检验公式根据以上公式,利用方差理论得到:如果x与y有线性关系,则其中F(1,n-2)表示第一自由度(分子的自由度)为1，第二自由度（分母的自由度）为n－2的F分布。105回归方程的显著性检验(F检验)F检验表(见书后附录P408)的使用若F<F0.05(1,n-2)则称x与y没有明显线性关系若F0.05(1,n-2)<F<F0.01(1,n-2)则称x与y有显著线性关系若F>F0.01(1,n-2)则称x与y有十分显著线性关系

106回归方程的显著性检验

(相关系数法)可作为考试题目利用S残及

证明|r|<=1从而由S总=S回+S残及S总=Lyy

S回=bLxy得到S残=Lyy-bLxy

因此|r|<=1107令:

于是

y=α+βx+ε可以写成矩阵形式:y=Xβ+ε且假定条件写成:E(ε)=0,D(ε)=σ2In回归模型的矩阵表示108回归模型的矩阵表示利用矩阵方法求回归系数β的最小二乘估计b

使得Q(b)=minQ(β)其中Q(β)=(y－Xβ)ˊ(y－Xβ)是1×1的矩阵，是多元函数109回归模型的矩阵表示补充知识:多元函数对向量求导规则:若多元函数y=f(x1,...xn)可以看作y为向量x=(x1,x2,..xn)‘的函数,记为y=f(x),且y对x的每个分量的微商都存在,则令称为y对向量x的微商.110回归模型的矩阵表示运算规则1.若y=x＇x

则2.若y=x＇Ax则

(推导从略)

特别地:若A为对称矩阵则:3.若y=a＇x111回归模型的矩阵表示求β的极值:112回归模型的矩阵表示若b是β的最小二乘估计,则有(X＇X)b=X＇y此式称为正规方程.如果X＇X可逆则可以得到b的解b=(X＇X)-1X＇y

113回归模型的矩阵表示验证b确实可以使Q达到最小值Q(β)=(y－Xβ)＇(y－Xβ)=(y－Xb+Xb－Xβ)＇(y－Xb+Xb－Xβ)=(y－Xb)＇(y－Xb)+(b－β)＇

X＇X(b－β)＋(y－Xb)＇X(b－β)

＋(b－β)＇

X＇(y－Xb)=(y－Xb)＇(y－Xb)+(b－β)＇

X＇X(b－β)＋(y－Xb)＇X(b－β)＋(b－β)＇

X＇(y－X(X＇X)-1X＇y)=(y－Xb)＇(y－Xb)+(b－β)＇

X＇X(b－β)＋(y－Xb)＇X(b－β)＋(b－β)＇[X＇y－X＇X(X＇X)-1X＇y]=(y－Xb)＇(y－Xb)+(b－β)＇

X＇X(b－β)＋(y－Xb)＇X(b－β)＋0=(y－Xb)＇(y－Xb)+(b－β)＇

X＇X(b－β)＋0＋0=Q(b)+(b－β)＇

X＇X(b－β)显然:

Q(β)>=Q(b),当且仅当β=b时Q(β)达到极小值.114回归模型的矩阵表示有关的性质证明：b是β的无偏估计115回归模型的矩阵表示相关性质证明：由正规方程，有b=(X'X)-1X'y

再根据协方差阵的性质：V(Ax)＝AV(x)A＇回归模型的矩阵表示Var（a）和Var（b）表达式a，b均是标量，不是向量117回归模型的矩阵表示Q(b)的表示：

其中是一个投影矩阵

(若A’=A,A2=A则A为投影矩阵)118残差分析所谓残差是指实际观察值与回归估计值的差，即

显然，有多少对数据，就有多少个残差。残差分析就是通过残差所提供的信息，分析数据的可靠性、周期性或其它干扰的一种数理统计方法。119残差分析

——残差的统计性质残差的期望为0由得到又由最小二乘估计中以及得所以120残差分析

——残差的统计性质《应用回归分析》张小蒂p94估计值与观测值的协方差残差分析

——残差的统计性质证明估计值的方差先证明：122残差分析

——残差的统计性质《残差分析

——残差的统计性质《再证：注意到：yi之间是独立的残差分析

——残差的统计性质《125残差分析

——残差的统计性质残差的方差：残差分析

——残差的统计性质《残差的方差中，σ一般是未知的，所以需要对σ进行估计。证明：E(Q(b))=(n-2)σ2S残=S总-S回S回=b2Lxx见PPT125先证：128129130考察残差图

张小蒂《应用回归分析》p84考察残差图

根据Guass-Markov条件

，回归模型中的随机扰动误差项εi~N(0,σ2)，如果模型对样本数据拟合良好的话，那么观测到的残差ei就应该反映εi的上述分布特性，即有ei~N(0,σ2)或者ei*=ei/σ~N(0,1)称为标准化残差。又由E(Q(b))=(n-2)σ2,可以得到标准化残差的估计量，且有P(|ei*|<2)=0.9545131考察残差图回归拟合良好：绝大多数数据在（-2，+2）水平区间，无任何系统趋势+2-20ei*

x

132回归函数具有曲线形式+2-20ei*

x

133样本数据中有异常点+20-2ei*

x

异常点134回归方程拟合不充分,较多数据在区间外,可能是回归模型的函数形式(对数,指数)选择不当,也可能是漏掉了重要自变量

+2-20ei*

x

135异方差,数据点往往出现系统变动趋势(下图是其中一种情形)可采用加权最小二乘法回归+2-20ei*

x

136误差项相关,也称自行关或序列相关,当数据为时间序列资料时,误差项自相关情况较多(这里只列举正相关情形εi>0的数据点对应εj>0的数据点)εiεj137残差分析

——预测预报根据历史数据回归得到方程后就可以用来做预报：假如根据新的数据x0代入回归方程得到新的预报值则有并且由于随机误差ε服从正态分布，残差也服从正态分布138残差分析根据分布理论可以得到：这里表示自由度为n-2的分布。并且有

139残差分析我们可以求得t1,t2,使覆盖y0的概率(t1<y0<t2)为1-α。考虑到，通常我们就取关于的对称区间。令140残差分析因此有：即141由可知：给定的αn越大，Lxx越大，x0越靠近x的均值，则Δ越小，预测精度越高。142残差分析实例：给定新值得所以，当已知当年积雪的最大深度为27.5尺时，以95%的概率断言灌溉面积在48.306千亩与56.12千亩之间。

143多元线性回归多元线性回归模型表示为：其中，是未知参数，而是m个可以精确测量并可控制的一般变量，是随机误差。和一元线性回归一样，常假定144多元线性回归为了建立回归方程，估计回归系数我们进行n次观察,得n组观察数据它们应有回归关系,可写成如下形式145多元线性回归设：是n组观察值是的估计量，记推导得146多元线性回归设：是β

的最小二乘解，则称为为正规方程用元素表示为：147多元线性回归对于正规方程如果满秩，即有逆矩阵存在。得最小二乘估计为

通常为了计算方便，往往并不先逆矩阵求b，而是通过解线性方程组(正规方程)来求b。148多元线性回归方程未知数较少时可用克莱姆法则求解在未知数较多时，可用高斯消去法等方法解方程149多元回归最小二乘估计的性质

——方阵“迹”的性质若A为n阶方阵，则方阵A的迹为其中λi为方阵A的特征根，或者其中aii

为方阵A的主对角线元素。（证明略）方阵的迹的基本性质：1 方阵正交变换后，迹不变 其中Q为正交阵2k为常数3若C为n×m阶矩阵，D为m×n阶矩阵，则有45150多元回归最小二乘估计的性质

现代地理统计分析证明：σ2的无偏估计为即证明：证明：已知151多元回归最小二乘估计的性质152多元回归最小二乘估计的性质容易得到：多元回归最小二乘估计的性质令：则有：154多元回归最小二乘估计的性质因此有：最后有：由此得的无偏估计为155多元回归最小二乘估计的性质(1/2)1．b是β无偏估计，这是因为

2.b的协方差矩阵为156多元回归最小二乘估计的性质对于线性回归方程系数的估计b，它的方差越小越好（方差越小表示波动越小，即Var（b）越小越好）在所有的线性无偏估计中最小二乘法得到的系数估计值方差最小无偏：

如果d为β的无偏估计，协方差阵记为Dd，若Dd-Db>=0(非负定)则称d的协方差阵Dd大于b的协方差阵Db

线性：对于β的一切线性组合c’β，若有Var(c’d)＞＝Var(c’b)

则称d（β的任一无偏估计）的协方差阵大于等于b（最小二乘估计）的协方差阵（c是常数向量）注意到c＇b=c'(XX)-1X'y

是y的线性函数（因此，b是线性估计量），因此可构造y的线性函数d'y作为c'β的任意一个线性无偏估计157多元回归最小二乘估计的性质3.对于线性回归模型，在c′β的一切线性无偏估计类中，c′b有最小方差（Guass-Markov定理）。证明:设d′y是c′β的任一线性无偏估计，则根据无偏估计的定义有因此有：对一切β成立，故必有这样由yi，yj之间的独立性（i≠j）158多元回归最小二乘估计的性质从而最后一步是因为I-P为投影阵，必为非负定阵。证毕。Guass-Markov定理说明：最小二乘估计在一切线性无偏估计中是最优的，因此也被称为最佳线性无偏估计159多元回归最小二乘估计的性质最小二乘估计存在的问题：β的估计的均方误差MSE(MeanSquaredErrors)160多元回归最小二乘估计的性质由于A>0（正定阵），因此A的所有特征根λi均为正数，且A-1的特征根为1/λi161多元回归最小二乘估计的性质由看出：当A的特征根中只要有一个根接近0，MSE会变大，回归系数估计值会变得很不稳定A的特征根接近0就是A接近奇异阵，导致X自变量存在线性关系当资料阵X的列向量之间存在近似线性相关关系时，称为复共线性关系（自变量不独立）对LS的改进：压缩估计（1960）、主成分估计（1965）、岭回归估计（1970）、特征根估计（1974）这些改进方法仍然是线性估计，但是不能保证无偏性11线性方程组病态的几何解释方程组是由多个线性方程组成每个线性方程都代表一条直线方程有解就是直线之间相交如果直线之间互相平行，那么或者有无数个交点（无穷多个解），或者没有交点（无解）直线之间如果接近平行，虽然有交点（有唯一解），但是直线参数略有小的变动都会导致交点大幅度变化，也就是说解是不稳定的，接近病态。162时间序列分析164时间序列的定义时间序列就是按照时间顺序排列的，随时间变化而且相互关联的数据序列。(如:年平均气温、降水量、海平面高程等)由于受到偶然因素的影响，每个时刻的取值是随机的不确定的，表现出某种随机性。165时间序列的定义能否把时间序列看成自变量是时间的一元回归分析？（不能直接使用，往往回归结果是一条水平线，相关系数为0没有意义，实际上我们用“自回归”研究时间序列）166时间序列的特点是一个随机变量的历史记录。自变量不限于时间，也可以是其他变量，但是必须是遵循某种顺序进行排列的数据的取值依赖于时间的变化，但是由于随机干扰的影响，不一定是时间t的严格函数(如同一元线性回归中xy没有严格函数关系)某一时刻上的取值具有一定的随机性，不可能完全准确地用历史值预测，某时刻在得到样本之前只是知道它的概率分布，而不知道具体的值前后时刻（不一定相邻）的数值具有一定的相关性时间序列往往呈现出一种趋势性（GDP）或周期性变化（海平面、气温）。167时间序列定义

——随机过程随机过程（stochasticprocess）的定义：定义：（从时间角度考察）若对于每个特定的t∈T（T是无穷集合，称为参数集）都存在一个随机变量X(t)，则称这一族的无穷多个随机变量是一个随机过程，表示为

{X(t),t∈T}或{Xt,t∈T}。对于某一固定的t，X(t)是一个普通随机变量，X(t)的取值空间称为相空间或状态空间。离散参数的随机过程也称为随机序列或时间序列。如果状态空间是离散的集合，参数t空间也是离散集合，则称此随机过程为链

时间tt状态空间t+kX（t）168随机过程与随机变量随机变量：描述随机现象 某班一天的出勤人数，某学院一天的耗电量

随机过程：描述随机现象的随时间的动态变化过程 某班每学期每天的出勤人数的情况，某学院每学期每天的耗电量 随机过程的现实（样本函数） 某班某学期每天的出勤人数的情况，某学院某学期每天的耗电量随机变量是随机过程的一个特例（固定时间点）。169随机过程的现实（样本函数）

Realization

170随机过程的统计特征

1.均值函数：给定随机过程{X(t),t∈T}，对于任意t∈T，若E[X(t)]存在，则称μt＝E[X(t)]，t∈T为随机过程{X(t),t∈T}的均值函数。均值函数是随机过程X(t)在各个时刻的摆动中心。2.方差函数

Var(Xt)=E[(Xt－μt)2]3.自协方差函数给定随机过程{X(t),t∈T}，取定t∈T，s∈T定义其自协方差函数为

γt,s=cov(Xt,Xs)=E[(Xt－μt)(Xs－μs)]当t=s时γt,t=Var(Xt)＝σt24.自相关函数：

ρt,s=γt,s

／(σtσs)171随机过程的统计特征时间tXs方差函数均值函数Xtμsμt172平稳随机过程

StationaryProcess平稳随机过程

n维分布函数F（当然所有统计特性）不随时间而变化的过程称为平稳随机过程（严平稳随机过程StrictlyStationaryProcess，狭义平稳随机过程） 为验证平稳性，需要计算所有有限维分布，通常在实际中难以满足以上要求。 实际上，通常平稳随机过程（序列）是指宽平稳随机过程（序列），定义如下：

173宽（弱）平稳过程

(序列)

WeaklyStationaryProcess(TimeSeries)设随机序列{X(t),t＝…,-2,-1,0,1,2,…}满足1.E[X(t)]＝μ＝常数；2.自协方差γt+k,t=γk=σk2与t无关(k＝…,-2,-1,0,1,2,…)则称Xt为（宽，弱）平稳随机序列（平稳时间序列），简称（宽，弱）平稳序列k实际上就是时间差，就是说相关函数与时间起点t无关，只与时间差值k有关。若当k=0时有γk＝σ2

，k≠0时γk＝0，则称Xt为平稳白噪声序列(WhiteNoise)174宽平稳过程与严平稳过程1）宽平稳随机过程不一定是严平稳随机过程。2）严平稳随机过程不一定是宽平稳随机过程如果自协方差存在则严平稳随机过程一定是宽平稳随机过程。

宽平稳要求自协方差和期望存在，而严平稳要求概率分布存在，并不断言期望和自协方差存在。比如柯西分布期望和方差均不存在。3）若随机过程是严平稳随机过程并且期望和自协方差都存在则一定是宽平稳，反之不一定成立4）对于正态随机过程（X(t)是正态分布），则严平稳与宽平稳是等价的

正态分布的期望和自协方差都存在，如果是严平稳则一定时宽平稳由于正态过程的概率密度由期望和自协方差完全决定的，因此，如果期望和自协方差不随时间平移变化，概率密度也不随时间平移变化，因此一个宽平稳的正态过程必然是严平稳的。175平稳序列的遍历性与参数估计1.如何知道一个时间序列是平稳序列？

Daniel检验 （Spearman相关系数）可参见<随机过程简明教程同济大学出版社何迎晖钱伟民51.716.53>P2002.对一个平稳序列，如何得到它的均值与自协方差函数和自相关函数?

均值求法：根据定义，可固定某个时间t，求各个现实在时刻t的平均。但是实际运用中存在困难。因为通常我们只能观察到一个现实。176平稳序列的遍历性与参数估计

3.平稳序列的遍历性对遍历性的直观解释：只要观测的时间足够长，随机过程的每个现实（样本函数）将“遍历”(统计参数以概率1相等)状态集中所有的取值情况。实际上，如果平稳序列存在遍历性，可以通过一个实现的样本数据就可以求出均值与协方差函数及自相关函数。因此，对多个现实的考察转换为对一个现实的考察，方便了统计特征的计算。177平稳序列的遍历性与参数估计当n充分大时，对于序列现实（样本）x1,x2,…,xn有γt,s=cov(Xt,Xs)=E[(Xt－μt)(Xs－μs)]μt＝E[X(t)]，t∈Tρt,s=γt,s

／(σtσs)均值函数自协方差函数自相关函数178平稳序列的遍历性与参数估计4.如何知道平稳序列具有遍历性?定理：如果平稳过程{X(t),t∈T}的相关函数满足，则{X(t),t∈T}具有均值各态遍历性。定理：如果平稳过程{X(t),t∈T}是具有0均值的正态过程，如果满足，则{X(t),t∈T}具有相关函数的各态遍历性工程应用中的实际做法是：先假定平稳过程具有各态遍历性，然后由此假定出发，对各种数据进行分析，在实践中考察是否会产生较大偏差，如果偏差较大，便认为此平稳过程没有各态遍历性。5.实际应用中哪些时间序列是平稳的并且具有遍历性的?大多数实际问题以及即将介绍的ARMA序列是具有遍历性的平稳序列。179平稳时间序列模型

若：Xt为零均值的平稳序列Xt－

Xt－1为一阶差分记做ΔXt，at

为随机扰动是平稳白噪声序列一阶自回归模型：

AR(1)：Xt＝φ1Xt－1+at

当φ1＝1时，称为随机游动，是AR(1)的特例表示为ΔXt=at说明系统的差异完全是由于外界的扰动造成的。n阶自回归模型AR(n)AutoRegressivemodel：

Xt＝φ1Xt－1+φ2Xt－2+φ3Xt－3+…+φnXt－n+atm阶移动平均模型MA(m)MovingAveragemodelXt＝at－θ1at－1－θ2at－2－θ3at－3－…－θmat－mn阶自回归m阶移动平均模型ARMA(n,m)AutoRegressiveMovingAveragemodelXt＝φ1Xt－1+φ2Xt－2+φ3Xt－3+…+φnXt－n+at－θ1at－1－θ2at－2－θ3at－3－…－θmat－m

180平稳时间序列模型平稳白噪声序列181上节课内容时间序列的定义随机过程的定义及相关概念随机过程与时间序列的关系随机过程的统计参数随机过程的各态遍历性及其检验随机过程的平稳性概念（重点掌握严平稳与宽平稳的关系）及其检验平稳时间序列常见模型（ARMAARMA）182本节课的主要内容求解AR（1）模型格林函数及其意义183为单摆建模

——平稳时间序列线性模型举例设第t个摆动周期的最大摆幅为Xt第t+1个摆动周期的最大摆幅为Xt+1考虑到阻尼系数ρ则有Xt+1=ρ

Xt又考虑到外界环境的其他随机影响at因此可建立模型：Xt+1=ρ

Xt+at其中{at}为白噪声序列，

|ρ|<1,t=…,-1,0,+1,…是一个一阶自回归模型AR（1)也称为Markov过程第t个摆动周期最大摆幅Xt184后移算子后移算子定义为：But=ut-1有如下性质：对于与时间t无关的随机变量u

，Bu=uB0≡1，称B0为恒等算子若c为一常数，则B(cut)=cBut对于任意两个序列ut和vt有B(ut±vt)=But±BvtBnut=ut-n后移算子在时间序列中的应用

求解AR(1)模型Xt＝φ1Xt－1+at首先引入后移算子B(Back)：

BXt＝Xt－1,BXt-1＝Xt－2,… Bat＝at－1,Bat-1＝at－2,…则AR(1)表示为Xt＝φ1BXt+at解之：显然代入差分方程可验证此式是原始差分方程的解变形为：得证。186格林函数对于一阶自回归模型，Gj＝φ1j说明AR(1)可用一个限阶的MA模型逼近。Gj是前j个时间单位之前进入到系统的扰动at－j

对现在行为的影响的权数。它反映了系统对干扰的响应的衰减的快慢程度。φ1越小衰减速度越快，完全由φ1所决定。方程系数函数φ1j

反映了系统对扰动的记忆程度，所以φ1j被称为记忆函数，又称格林(Green’sfunction)一般用Gj表示。187格林函数与平稳性AR(1)系统的平稳性|φ1|<1平稳 随着时间的推移扰动项的影响以指数方式衰减，系统最终趋向平衡位置|φ1|＝1临界平稳 扰动项的影响没有衰减，但由于扰动是白噪声，系统的状态是有界的|φ1|>1非平稳 扰动项的影响在放大，系统不断以指数方式加速远离平衡位置，且永远不会恢复到平衡位置平衡位置188Xt方差

|φ1|＝1时方差不存在，不是宽平稳，因此随机游动Xt＝

Xt－1+at不是平稳的 随机游动的平稳性189Wold分解

——从线性空间的角度解释格林函数

由于随机扰动at是相互独立的，所以可以看作线性空间的一组基（相互之间没有线性关系），格林函数就是关于基的坐标。Xt就是这个空间的一个点。因此格林函数Gi也称Wold系数，

也叫做Wold分解式。190时间序列时域分析工具

——线性常系数差分方程常系数线性微分方程描述连续时间系统动态性，常系数线性差分方程描述离散时间系统的动态性。线性常系数差分方程的解法可比拟微分方程的解法，实际上，只不过差分方程是离散的，微分方程是连续的。191线性常系数差分方程n阶差分方程：y(k+n)+an-1y(k+n-1)+an-2y(k+n-2)+…+a0y(k)=u(k)其中，ai为系统参数的函数，当ai为常数时，就是常系数n阶差分方程。u(k)是个离散序列，也称作驱动函数，y(k+n)是系统的响应。当u(k)＝0时，称方程y(k+n)+an-1y(k+n-1)+an-2y(k+n-2)+…+a0y(k)=0为齐次常系数差分方程。192线性常系数差分方程的解求解线性常系数差分方程就是在给定n个初始条件y(0)，y(1)，y(2)，…，y(n-1)求出y(n),y(n+1)…对于n阶自回归m阶移动平均模型ARMA(n,m)

Xt＝φ1Xt－1+φ2Xt－2+φ3Xt－3+…+φnXt－n+at－θ1at－1－θ2at－2－θ3at－3－…－θmat－m就是给定Xt－1

，Xt－2

，Xt－3

，…，

Xt－n

求出Xt其中at－θ1at－1－θ2at－2－θ3at－3－…－θmat－m=u(t)可以看作驱动函数求解差分方程与求解微分方程类似：先求出对应齐次方程的通解，然后求出原方程的一个特解，特解和通解的线性组合就构成了原方程的解193线性常系数差分方程及其解的一般形式求解步骤：1.求出相应齐次差分方程的通解设：Y(k)=λk

是齐次差分方程y(k+n)+an-1y(k+n-1)+an-2y(k+n-2)+…+a0y(k)=0的一个解，则必有：λk+n+an-1λk+n-1+…+a0λk=0;从而有特征方程：λn+an-1λn-1+…+a0=0;求得方程的n个特征根λi(i=1,..,n)即可得到齐次差分方程的通解Y(k)=ΣCiλik其中，Ci为任意实数，λi可能是实数或复数。2.求一个原方程的特解一般令y(k)=i（i是常数）即可求得特解。3.原方程的解就是通解与特解的线性组合194线性常系数差分方程求解举例求解二阶非齐次差分方程解：求出对应齐次方程的通解。设是对应齐次方程的一个解，则有：解之得：因此得到通解：195线性常系数差分方程求解举例

求出原方程的特解：令：代入原方程：得到：特解为：原方程的解：196格林函数的解

——隐式解ARMA(2,1)模型：Xt－φ1Xt－1－φ2Xt－2＝at－θ1at-1的解设为则用B算子表示为代入模型方程得到：197

若等式成立则对应同次幂的系数必然相等于是得到：在系统参数已知情况下，根据以上表达式就可以递推计算出所有的Gj

，当j充分大时即可计算出方程的解。198格林函数的解

——显式解ARMA(2,1)是一个二阶非齐次差分方程：Xt－φ1Xt－1－φ2Xt－2＝at－θ1at-1解之：求通解：将Xt－2＝λk代入上式得到：λk+2

－φ1λk+1－φ2λk=0即得特征方程：λ2

－φ1λ－φ2=0求解特征方程得到通解：Gj＝g1λ1j

＋g2λ2j其中g1，g2是任意常数。199格林函数显式解求解系数g1，g2根据隐式解给出的初始条件：G0＝1；G1－φ1G0＝－θ1；有：根据韦达定理有：因此有：解之得：200格林函数显式解因此得到ARMA(2,1)系统的格林函数为：存在两个共轭复数根的情况（略）201格林函数ARMA(1,2)是一个二阶非齐次差分方程：at－θ1at-1－θ2at-2＝Xt－φ1Xt－1解之：求通解：将at－2＝vk代入上式得到：vk+2

－θ1vk+1－θ2vk=0即得特征方程：

v2

－θ1v－θ2=0求解特征方程得到通解：Ij＝g1v1­j

＋g2v2­j其中g1，g2是任意常数。202格林函数因此得到ARMA(1,2)系统的格林函数为：与格林函数相似,只有当|v1|<1且|v2|<1系统才是可逆的(Ij→0)。存在两个共轭复数根的情况（略）203格林函数的解

——B算子河海P171）考虑模型AR(2)用后移算子表示为令则有考虑到格林函数Gj得到格林函数由上式可知当j→∞时有Gj→0，因此该系统是平稳的

205求解格林函数-B算子

（P59）考虑模型ARMA(2,1)用后移算子表示为则有

可见：系统要稳定(级数收敛)要求|λ|<1,因此系统稳定的条件可描述为|1/λ|>1,而1/λ是φ(B)=0的根，或者说要求B的所有根要在单位圆外部207得到格林函数Gj又由于得到：即可解得λ，可见λ只是φ的函数求解格林函数

（P59）208格林函数ARMA(2,1)系统的稳定性：由知只有当j→∞时有Gj→0，系统是平稳的。由格林函数Gi＝g1λ1j

＋g2λ2j

知只有当|λ1|<1且|λ2|<1系统才是平稳的。因此：由于系统的平稳性只与自回归参数有关，与移动平均参数无关。（可从格林函数求解的表达式看出来）因此，所有ARMA(2,m)系统的平稳性条件都相同，都是上式。209格林函数-1-0.500.51-2-1012φ1φ2可以得到系统的平稳区域：210逆函数和可逆性定义：如果一个过程可以用一个无限阶的自回归模型逼近，即逆函数存在，就称该过程具有可逆性。AR(n)的逆函数（略）MA(1)的逆函数：因此：可逆性条件：才能保证211

格林函数与逆函数对于模型令：则原式变为为：因此，逆函数与格林函数之间符号相反，参数互换即可得到对方。（具有对偶性）212格林函数与逆函数1.格林函数的平稳性只与AR模型的特征根λk有关;2.逆函数的可逆性仅与MA模型的特征根vk有关.3.对于ARMA模型可逆又平稳的条件是特征方程的所有根满足:|λk|<1且|vk|<14.判断平稳性与可逆性使用的方法是完全一致的,只是变量记号不同对于高阶的ARMA通过解方程的方法求解逆函数和格林函数是困难的,实