自然正交函数分析(EOF)程序
EOF分析及其应用
EOF分析及其应用
一、EOF分析是什么
EOF分析（Empirical Orthogonal Function Analysis）是一种常用
的时间-空间统计分析方法，它是由把空间上的一维观测或多维观测数据
矩阵投影到一个更特别的模型空间中，然后对该模型空间中的变换数据进
行分析从而推算出有关的特征参数的一种分析方法。
二、EOF分析的原理
EOF分析由英国天文学家Harold E. Jeffreys （1891-1989）于
1931年提出。
它利用最小二乘估计法，把空间上一维或多维观测的数据
矩阵投影在一个特定的模型空间中，然后对该模型空间中变换的数据进行
分析，从而推算出有关的特征参数。
EOF分析的核心理论是“变换空间”，即给定一个多维空间Vn，找出一个低维变换空间Vm具有一定的特殊性质（如基Vm上的每一列向量的模具有最小值，它们张成一个最小的模型空
间上），使得数据在其中具有最好的表示，且在该变换空间中可以表示出
空间统计分布的特性。
三、EOF分析的应用
（1）短时间强对流预报
短时间强对流预报是一种有效的大气环境监测技术，它依据大气各层
能量释放特征进行短时间的天气预报。
EOF方法运用了空间观测数据，可
以对大气能量释放做出准确的模拟分析，从而预测出未来几小时内这一区
域内的强对流天气预报。
（2）大气环流异常研究。
第五章主成分分析（2）（主成分回归经验正交分解EOF）
5.4 主成分聚类与主成分回归5.4.1 变量聚类与样品分类主成分分析可用于聚类：变量聚类与样品聚类。
变量聚类：由主成分系数的差异，可将变量聚类。
例如例5.5中第2主成分中murder,rape, assult系数为负的， burglary,larceny, auto系数是正的。
按系数正负可把7个变量分为两类： murder, rape, assult属于暴力程度严重的一类；burglary,larceny,auto属于暴力程度较轻的一类。
按照这种方法，根据主成分系数的正负可以将变量聚类。
样品聚类：如果2个主成分能很好的概括随机向量的信息，计算每个样品的这两个主成分得分，把他们的散点图画出来，就能从图上将样品分类。
例5．5（续2）按照第一、第二主成分得分，画出散点图data crime; /*建立数据集crime*/input state $ 1-15 murder rape robbery assult burglary larceny auto;/*建立变量state murder rape robbery assult burglary larceny auto。
state $ 1-15表示前15列存州名。
murder rape robbery assult burglary larceny auto 表7种罪的犯罪率*/cards; /*以下为数据体*/Albama 14.2 25.2 96.8 278.3 1135.5 1881.9 280.7Alaska 10.8 51.6 96.8 284.0 1331.7 3369.8 753.3Arirona 9.5 34.2 138.2 312.3 2346.1 4467.4 439.5Arkansas 8.8 34.2 138.2 312.3 2346.1 4467.4 439.5Califonia 11.5 49.4 287.0 358.0 2139.4 3499.8 663.5Colorado 6.3 42.0 170.7 292.9 1935.2 3903.2 477.1Conecticat 4.2 16.8 129.5 131.8 1346.0 2620.7 593.2Delaware 6.0 24.9 157.0 194.2 1682.6 3678.4 467.0Florida 10.2 39.6 187.9 449.1 1859.9 3840.5 351.4Geogia 11.7 31.1 140.5 256.5 1351.1 2170.2 297.9Hawaii 7.2 25.5 128.0 64.1 1911.5 3920.4 489.4Idaho 5.5 19.4 39.6 172.5 1050.8 2599.6 237.6Illinois 9.9 21.8 211.3 209.0 1085.0 2828.5 528.6Indiana 7.4 26.5 123.2 153.5 1086.2 2498.7 377.4Iowa 2.3 10.6 41.2 89.8 812.5 2685.1 219.9Kansas 6.6 22.0 100.7 180.5 1270.4 2739.3 244.3Kentaky 10.1 19.1 81.1 123.3 872.2 1662.1 245.4Loisana 15.5 30.9 142.9 335.5 1165.5 2469.9 337.7Maine 2.4 13.5 38.7 170.0 1253.1 2350.7 246.9Maryland 8.0 34.8 292.1 358.9 1400.0 3177.7 428.5Masschusetts 3.1 20.8 169.1 231.6 1532.2 2311.3 1140.1Michigan 9.3 38.9 261.9 274.6 1522.7 3159.0 545.5Minnesota 2.7 19.5 85.9 85.8 1134.7 2559.3 343.1Mississippi 14.3 19.6 65.7 189.1 915.6 1239.9 144.4Missouri 9.6 28.3 189.0 233.5 1318.3 2424.2 378.4Montana 5.4 16.7 39.2 156.8 804.9 2773.2 309.3Nebraska 3.9 18.1 64.7 112.7 760.0 2316.1 249.1Nevada 15.8 49.1 323.1 355.0 2453.1 4212.6 559.2Mew Hampashare 3.2 10.7 23.2 76.0 1041.7 2343.9 293.4New Jersey 5.6 21.0 180.4 185.1 1435.8 2774.5 511.5New Maxico 8.8 39.1 109.6 343.4 1418.7 3008.6 259.5New York 10.7 29.4 472.6 319.1 1728.0 2782.0 745.8North Carolina 10.6 17.0 61.3 318.3 1154.1 2037.8 192.1North Dakoda 100.9 9.0 13.3 43.8 446.1 1843.0 144.7Ohio 7.8 27.3 190.5 181.1 1216.0 2696.8 400.4Oklahoma 8.6 29.2 73.8 205.0 1288.2 2228.1 326.8Oregan 4.9 39.9 124.1 286.9 1636.4 3506.1 388.9Pennsyvania 5.6 19.0 130.3 128.0 877.5 1624.1 333.2Rhode Island 3.6 10.5 86.5 201.0 1849.5 2844.1 791.4South Carolina 11.9 33.0 105.9 485.3 1613.6 2342.4 245.1South Dakoda 2.0 13.5 17.9 155.7 570.5 1704.4 147.5Tennessee 10.1 29.7 145.8 203.9 1259.7 1776.5 314.0Texas 13.3 33.8 152.4 208.2 1603.1 2988.7 397.6Utah 3.5 20.3 68.8 147.3 1171.6 3004.6 334.5Vermont 1.4 15.9 30.8 101.2 1348.2 2201.0 265.2Virginia 9.0 23.3 92.1 165.7 986.2 2521.2 226.7Wasinton 4.3 39.6 106.2 224.8 1605.6 3386.9 360.3West Viginia 6.0 13.2 42.2 90.9 597.4 1341.7 163.3Wiskonsin 2.8 12.9 52.2 63.7 846.9 2614.2 220.7Wyoming 5.4 21.9 39.7 173.9 811.6 2772.2 282.0;proc princomp out=crimprin n=2;var murder rape robbery assult burglary larceny auto;run;PROC PLOT data=crimprin;PLOT PRIN2*PRIN1=STATE/VPOS=31;TITLE2 ‘PLOT OF THE FIRST TWO PRINCIPAL COMPONENTS’;RUN;例5．7 （气温分析）本例的输入资料文件（TEMPERA T）是美国六十四个城市一月与七月的平均日温。
常用数据分析方法介绍
样本长度、时间尺度个数、起始时间 尺度、时间尺度间距
参数说明
• （4）小波分析程序输出结
年份
时间尺度
果文件为WA文件夹下的
“Fileout.txt”，给出了年份
小波系数
、时间尺度以及小波系数
值；
20
18
16
时 14 间 尺 12 度 10 /8 a6
4
2 1961 1964 1967 1970 1973 1976 1979 1982 1985 1988 1991 1994 1997 2000 2003 2006 2009
天长
界首 临泉
太和 阜阳
阜南
涡阳
宿县
灵壁 泗县
蒙城 利辛
固镇
五河
颍上
怀远蚌埠 凤阳
凤台 淮南 寿县 长丰
定远
明光
霍邱
来安 滁州
天长
全椒
金寨
六安 霍山
岳西
合肥 肥东
肥西 舒城
含山 和县马鞍山
巢湖
当涂
庐江 桐城
无为
芜湖
铜陵
繁昌
芜湖县
南陵
宣城
郎溪 广德
潜山 太湖 怀宁
宿松
望江
枞阳 贵池
青阳
安庆
九华
泾县
东至
Fortran计算程序中需要修改的参数主要有：N（样本长度）、 NYEAR（起始年份）
样本长度、起始年份
• （4）MK检验程序输出结果文 件为MK文件夹下的 “Fileout.txt”，其中第一列为 年份；第二、三列分别为UF 和UB统计量值；第四、五列 为显著性水平。
年份
UF
UB
显著性水平
自然正交函数分析程序
自然正交函数分析程序EOF方法是基于假设，即数据可以被表示为一系列正交函数的线性组合。
这些正交函数称为EOF模态，并代表了数据中的主要模式。
每个EOF 模态都具有相应的权重，称为EOF系数，用于描述该模态在总方差中的贡献程度。
EOF方法的步骤如下：1.数据预处理：首先，要对原始数据进行预处理。
这可以包括去除重复数据、去除异常值、进行数据平滑处理等。
2.协方差矩阵计算：接下来，需要计算数据的协方差矩阵。
协方差矩阵描述了数据中不同维度之间的相关性。
3.特征值分解：通过对协方差矩阵进行特征值分解，可以得到特征值和特征向量。
特征值表示了每个特征向量对总方差的贡献程度。
4.选择模态：根据特征值的大小，可以选择保留最重要的EOF模态，从而降低数据维度。
5.计算EOF系数：对于每个选定的EOF模态，可以计算其相应的EOF 系数，用于描述该模态在数据中的贡献程度。
6.重构数据：最后，通过将所有选定的EOF模态与相应的EOF系数进行线性组合，可以重构原始数据。
这样可以去除一些噪音和次要特征，从而提取出原始数据中的关键模式。
EOF方法有许多应用，特别是在气候学、地球物理学和图像处理等领域。
在气候学中，EOF可以帮助我们理解地球上不同地区的温度、降水和风向等变化模式。
在地球物理学中，EOF可以帮助我们分析地震数据、地磁数据和重力场数据等。
在图像处理中，EOF可以帮助我们提取图像中的关键特征，用于图像分类和识别。
总之，自然正交函数分析（EOF）是一种强大的数学工具，用于处理具有时间或空间相关性的数据。
通过对数据进行正交分解，EOF可以提取出关键的时间或空间模式，并帮助我们理解和分析数据中的重要特征。
EOF应用：从数据预处理到详细分析
By lqouc
1. 什么是 EOF，它的作用是什么。 1.1 什么是 EOF 关于 EOF 要先从主成分分析说起，主成分分析是多元统计分析中重 要的一部分，是一种从多个变量化为少数变量的统计方法，利用多个 变量之间相互关系构造一些新的变量， 这些新的变量不仅能综合反映 原来多个变量的信息，而且彼此之间是相互独立的，同时是按方差贡 献大小排列的，这种统计处理方法称为主成分分析。主成分分析在气 象应用中称为经验正交函数（EOF）分解。 1.2EOF 的用途 对于一个气象要素， 我们通常有 m 个空间点或者台站， 有 n 次观测， 这样组成的矩阵中的任意元素就表示了某一空间某一时刻的函数， 我 们希望能将这样的时空函数分解成空间函数与时间函数两部分的线 性组合。根据主成分的性质，主成分是按其方差贡献大小排列的，而 且是相互独立的， 那么可以用前几个时间函数与对应的空间函数的线 性组合，对原始场做出估计和解释，这就是经验正交函数分解的主要 目的。 2. EOF 的数据预处理 EOF 只是个统计学的方法， 本身不带有任何物理意义， 更不会揣摩作 者的意图，所以在数据导入之前需要对数据进行分析和预处理。以免 得到错误的或者不理想的结果。 在此处所说的预处理不是指一般 EOF
小时间序列的自由度。3.带通滤波也是常用的方法（本人没用过） ， 其优点是可以选定一定的频率范围，缺点是边界处处理不是很清晰。 4.谐波滤波，以傅里叶函数为基函数对时间序列进行逼近，其优点是 可以较准确的得到选取的频段信号，缺点是选的基函数有局限性，而 且结果和时间序列的长度有关。5.线性去趋势可以去除时间序列的线 性趋势信号，但是需要这一线性趋势通过显著性检验。 2.3 如何合理选定分析对象 上面谈到的是滤波的方法， 但是如果我们的数据是一些大家不熟悉的 数据，我们并不知道它都主要包含何种尺度的信号，也不知道各个主 要尺度信号的强弱，那就需要先对时间序列进行分析。对于时间序列 的分析，我们可以采用 1.谐波滤波，看各个频率的数值大小。2.功率 谱分析，得到显著周期。3.小波分析，同样可以得到时间序列的多尺 度变化特征。 在此，我推荐的方法是结合空间利用方差分析，因为以上的分析我们 都是忽略了空间的影响， 一种要素的时间变化特征是会随着空间变化 的。例如，对中国地区做某一要素的 EOF 分析，得到的结果不能通 过检验（检验的方法，后面再说） ，这个时候我们就需要考虑是否一 些地区的目标信号不强，而另外一些地区目标信号很强，这样的话就 只需要分析目标信号很强的地区，即只对特定区域进行 EOF 分析。 结合空间的方差分析， 首先需要对要素每一个空间点的时间序列进行 滤波，得到各个不同频率的信号（从季节内到线性趋势） 。对每个平 率的信号求方差，得到了各个频率的方差的空间分布。在分析的过程
EOF应用：从数据预处理到详细分析
响因子， 进行简单相关、 复相关和偏相关分析， 确定可能的影响因子。 确定了影响因子之后可以尝试用多元回归分析， 探讨这些因子与研究 要素之间的可预报性。 除了以上提到的分析，还可以根据自己的目的增加分析的内容。 5. 不同类型的 EOF 5.1EOF 本身的变化 对于 EOF 的介绍很多的参考书籍都将其用于时空分离，也就是用在 了空间和时间构成的三维场。但是实际上，我们回归最前面的 EOF 的出处， 可以看出最原本的主成分分析并没有限定要素是时空的函数。 这种方法只是通过引入新变量来达到数组降维的效果。 所以我们可以 在应用中进行多种尝试，只要能在物理上找到合理的解释就没问题。 因为，这终究只是一种数学工具。 举个例子，我们将一个 30 年长度月分辨率的时间序列，写成一个 30*12 的数组，第一维 30 年，第二维是 12 个月，这样以 30 年为我 们通常认为的时间，12 个月为‘空间’ ，进行 EOF 分析，得到的结果 可以揭示不同模态下 12 个月分别在这 30 年中的变化。 除此之外还有很多种用法，在此不再赘述，仅作抛砖引玉。 5.2 多变量 EOF（MV-EOF） EOF 分析时， 不仅会研究某一要素的时空特征， 有时也会研究某现象 的时空特征，而这些现象往往不能用单一的要素来表征，这时候就需 要用到了多变量的 EOF。 例如，研究海洋大陆的季风系统时空变化特征，很可能要考虑到
小时间序列的自由度。3.带通滤波也是常用的方法（本人没用过） ， 其优点是可以选定一定的频率范围，缺点是边界处处理不是很清晰。 4.谐波滤波，以傅里叶函数为基函数对时间序列进行逼近，其优点是 可以较准确的得到选取的频段信号，缺点是选的基函数有局限性，而 且结果和时间序列的长度有关。5.线性去趋势可以去除时间序列的线 性趋势信号，但是需要这一线性趋势通过显著性检验。 2.3 如何合理选定分析对象 上面谈到的是滤波的方法， 但是如果我们的数据是一些大家不熟悉的 数据，我们并不知道它都主要包含何种尺度的信号，也不知道各个主 要尺度信号的强弱，那就需要先对时间序列进行分析。对于时间序列 的分析，我们可以采用 1.谐波滤波，看各个频率的数值大小。2.功率 谱分析，得到显著周期。3.小波分析，同样可以得到时间序列的多尺 度变化特征。 在此，我推荐的方法是结合空间利用方差分析，因为以上的分析我们 都是忽略了空间的影响， 一种要素的时间变化特征是会随着空间变化 的。例如，对中国地区做某一要素的 EOF 分析，得到的结果不能通 过检验（检验的方法，后面再说） ，这个时候我们就需要考虑是否一 些地区的目标信号不强，而另外一些地区目标信号很强，这样的话就 只需要分析目标信号很强的地区，即只对特定区域进行 EOF 分析。 结合空间的方差分析， 首先需要对要素每一个空间点的时间序列进行 滤波，得到各个不同频率的信号（从季节内到线性趋势） 。对每个平 率的信号求方差，得到了各个频率的方差的空间分布。在分析的过程
eof经验正交函数
eof经验正交函数正交设计是一种在实验设计中应用广泛的方法，它可以有效地降低实验的复杂性并提高实验结果的可信度。
在正交设计中，经验正交函数（EOF）是一种重要的工具，它可以帮助研究人员选择合适的实验因素和水平，以获得准确且可靠的实验结果。
经验正交函数是一组具有正交性质的基础函数，它们可以表示实验因素的不同水平。
通过使用经验正交函数，研究人员可以在有限的实验次数下，对多个因素进行全面的测试，从而节省时间和资源。
在实际应用中，经验正交函数可以用于设计和优化各种工程和科学实验。
例如，在材料科学领域，研究人员可以使用经验正交函数来优化材料的物理和化学性质。
在制造业中，经验正交函数可以用于优化生产过程中的参数设置，以提高产品质量和生产效率。
经验正交函数的选择和使用需要考虑因素的数量和水平，以及实验的目标和约束条件。
通常，研究人员可以使用经验正交函数表来选择合适的函数和水平组合。
经验正交函数表是经过统计和数学分析得出的，可以帮助研究人员快速准确地选择合适的函数。
除了选择适当的经验正交函数，研究人员还需要确定实验的因素和水平。
因素是指影响实验结果的变量，而水平则是指每个因素的具体取值。
通过选择合适的因素和水平，研究人员可以在实验设计中获得准确和可靠的结果。
在进行实验时，研究人员需要根据经验正交函数和选择的因素水平制定实验计划。
实验计划应考虑到实验次数、实验顺序和实验条件等因素，以确保实验结果的可靠性和可重复性。
经验正交函数的应用可以帮助研究人员在有限的实验条件下，获得准确和可靠的实验结果。
它不仅可以提高实验效率，还可以降低实验成本。
因此，经验正交函数在工程和科学领域中得到了广泛的应用和重视。
经验正交函数是一种重要的实验设计工具，可以帮助研究人员选择合适的实验因素和水平，以获得准确且可靠的实验结果。
它在各个领域的应用都取得了显著的成果，并为实验设计提供了有效的方法和策略。
通过合理地运用经验正交函数，研究人员可以在实验中更好地发现和理解因果关系，为工程和科学领域的发展做出贡献。
eof分析
事实上，这种想法是可以实现的，主分量 分析方法就是综合处理这种问题的一种强有力 的工具。 主分量分析是把原来多个变量划为少数几 个综合指标的一种统计分析方法。 从数学角度来看，这是一种降维处理技术。
在实际问题研究中，为了全面、系统地分析问 题，我们必须考虑众多影响因素。这些涉及的 因素一般称为指标，在多元统计分析中也称为 变量。因为每个变量都在不同程度上反映了所 研究问题的某些信息，并且指标之间彼此有一 定的相关性，因而所得的统计数据反映的信息 在一定程度上有重叠。在用统计方法研究多变 量问题时，变量太多会增加计算量和增加分析 问题的复杂性，人们希望在进行定量分析的过 程中，涉及的变量较少，得到的信息量较多。 主成分分析正是适应这一要求产生的，是解决 这类题的理想工具。
n
ki
xi )( x kj x j )
2
( xki xi )
( x kj x j ) 2
k 1
n
（1.3.2）
（二）计算特征值与特征向量
① 解特征方程 I R 0 ，常用雅可比法 （Jacobi）求出特征值，并使其按大小顺序排 列 1 2 p 0 ； ② 分别求出对应于特征值 i 的特征向量
主分量分析与核主分量分析
第一节 主分量分析
第二节 核主分量分析
第一节 主分量分析
概 述 主分量分析的基本原理 主分量分析的计算步骤 主分量分析主要的作用 主分量分析方法应用实例
一、概述
许多系统是多要素的复杂系统，多变量问 题是经常会遇到的。变量太多，无疑会增加分 析问题的难度与复杂性，而且在许多实际问题 中，多个变量之间是具有一定的相关关系的。 因此，人们会很自然地想到，能否在相关 分析的基础上，用较少的新变量代替原来较多 的旧变量，而且使这些较少的新变量尽可能多 地保留原来变量所反映的信息？
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
５．３自然正交函数分析（ＥＯＦ）程序近年来，自然正交函数（又称经验正交函数）展开在气象上应用比较广泛。
这种正交函数展开不象三角函数展开、球函数展开那样有固定的展开形式。
它无固定的函数形式，不是事先人为地给定典型场函数，图形是由场本身来决定的，它具有收敛快又能更好地反映出场的基本结构的特征。
它可以在有限的区域中进行，既可以取空间不同站点进行分解，也可以对同一站点的不同时间、不同高度的多种要素进行综和分析。
因此它在气象中具有广泛的应用，可用于气象要素场分析、大气垂直结构分析、动力模型垂直分层等。
5.3.1功能计算要素场的自然正交函数分解。
5.3.2方法说明自然正交函数分解是针对气象要素场进行的，它的基本思想是把包含ｐ个空间点（或ｐ个变量）的n 个时次的观测场随时间进行分解，即将某一区域的气象要素场序列ij F （i=1, 2, …,p ；j=1,2,…,n ，即p 个空间点的n 个时次的观测资料）分解成相互正交的时间函数与相互正交的空间函数的乘积之和，常把空间函数ik v 看作典型场，时间函数kj t 看作典型场的权重系数，则不同时间的要素场是若干个典型场按不同权重线性叠加的结果，各个场之间的差别就在于各典型场的系数不同。
则气象要素场可以表示为∑=+++==p 1k pj ip j 22i j 11i k j ik ij t v t v t v t vF （５．３．１）其中F ij 表示第i 个场中的第j 个测点的观测值。
可将(5.3.1)是写为矩阵的形式VT F = （５．３．２）式中F 为n p ⨯阶的均值为0的资料阵，V 为p p ⨯阶的空间函数阵，T 为n p ⨯阶的时间函数阵。
由于V 和T 是根据场的资料阵F 进行分解而得到的，分解的函数没有固定的函数形式，因而称为“经验”的，另外，我们还要求这种分解具有“正交”性，即要求满足下式⎪⎪⎭⎪⎪⎬⎫≠=='≠=='∑∑==)l k (0t t t t )l k (0v v v v n 1j lj k j l k p 1i il ik l k （5．3．３）事实上，我们对（５．３．２）式右乘T '可得V T VT F F ''=' （５．３．４）因F F '是p p ⨯阶对称阵，其元素为距平变量的交叉积。
根据实对称矩阵的分解定理有 V V ΛF F '=' （５．３．５）其中Λ是F F '矩阵的特征值组成的对角阵，V 是对应的特征向量为列向量组成的矩阵。
比较（５．３．４）和（５．３．５）式可知ΛT T =' （５．３．６）又根据特征向量的性质有I V V V V ='=' （５．３．７）式中为I单位矩阵。
显然（５．３．６）和（５．３．７）式满足（５．３．３）式的要求。
由F'矩阵的特征向量求得，而时间函数则可利用（５．３．２）式左此可知空间函数矩阵可从F乘V'得到，即T'=（５．３．８）VF5.3.3子程序语句CALL EOF(X,P,N,XF)5.3.4哑元说明X——输入变量，二维实型数组，大小为P⨯N，存放原始观测值。
P——输入整型变量，空间格点数。
N——输入整型变量，序列的时间长度。
XF——输出变量，二维实型数组，大小为P⨯N，存放恢复值。
5.3.5子程序SUBROUTINE EOF(X,P,N,LW,XF)INTEGER::PINTEGER::NINTEGER::LWREAL(8),DIMENSION(P,N)::X,XFREAL(8),DIMENSION(P,P)::A,V,V1REAL(8),DIMENSION(P,N)::TREAL(8),DIMENSION(P)::B,GM,GAREAL(8),DIMENSION(P,LW)::VFREAL(8),DIMENSION(LW,N)::TF! 求X乘以X的转置,即A=XXˊDO I=1,PDO J=1,PA(I,J)=0DO K=1,NA(I,J)=A(I,J)+X(I,K)*X(J,K)END DOEND DOEND DO! 用Jacobi法求A的特征值和特征向量! 返回时B存放矩阵的全部特征值，V存放特征向量为列组成的矩阵CALL JCB(A,P,1.0E-6,V,B,L)DO I=1,PGA(I)=0DO J=1,IGA(I)=GA(I)+B(J)END DOEND DODO I=1,PGM(I)=GA(I)/GA(P)END DODO I=1,PV1(I,J)=V(J,I)END DOEND DOT=MATMUL(V1,X)WRITE(12,'(" 特征值")')WRITE(12,'(<P>I10)')(I,I=1,P)WRITE(12,'(3X,<P>D10.4)')BWRITE(12,'(" 解释的方差(%)")')WRITE(12,'(<P>I7)')(I,I=1,P)WRITE(12,'(3X,<P>F7.2)')GM*100WRITE(12,'(" 特征向量为列组成的矩阵，即空间函数V")')WRITE(12,'(<P>F7.4)')((V(I,J),J=1,P),I=1,P)WRITE(12,'(" 时间函数T")')WRITE(12,'(<P>F10.4)')((T(I,J),J=1,P),I=1,N)DO I=1,PDO J=1,LWVF(I,J)=V(I,J)END DOEND DODO I=1,LWDO J=1,NTF(I,J)=T(I,J)END DOEND DOXF=MATMUL(VF,TF)ENDSUBROUTINE JCB(A,N,EPS,V,B,L)! A:调用时存放实对称矩阵! B:返回时存放矩阵的全部特征值! V:存放特征向量,其中第i列为与第i个特征值相对应的特征向量! EPS:存放精度要求REAL(8),DIMENSION(N,N)::A,VREAL(8),DIMENSION(N)::BREAL(8)::FM,CN,SN,OMEGA,X,YINTEGER::P,QL=1V=0DO I=1,NV(I,I)=1END DO10 FM=0DO I=1,NIF(I/=J.AND.ABS(A(I,J))>FM)THENFM=ABS(A(I,J))P=IQ=JEND IFEND DOEND DOIF(FM<EPS)THENL=1GOTO 15END IFIF(L>100)THENL=0WRITE(12,'("L=",I3,2X,"没有达到精度要求")')L GOTO 15END IFL=L+1X=-A(P,Q)Y=(A(Q,Q)-A(P,P))/2OMEGA=X/SQRT(X*X+Y*Y)IF(Y<0)OMEGA=-OMEGASN=1+SQRT(1-OMEGA*OMEGA)SN=OMEGA/SQRT(2*SN)CN=SQRT(1-SN*SN)FM=A(P,P)A(P,P)=FM*CN*CN+A(Q,Q)*SN*SN+A(P,Q)*OMEGA A(Q,Q)=FM*SN*SN+A(Q,Q)*CN*CN-A(P,Q)*OMEGA A(P,Q)=0A(Q,P)=0DO J=1,NIF(J/=P.AND.J/=Q)THENFM=A(P,J)A(P,J)=FM*CN+A(Q,J)*SNA(Q,J)=-FM*SN+A(Q,J)*CNEND IFEND DODO I=1,NIF(I/=P.AND.I/=Q)THENFM=A(I,P)A(I,P)=FM*CN+A(I,Q)*SNA(I,Q)=-FM*SN+A(I,Q)*CNEND IFEND DOFM=V(I,P)V(I,P)=FM*CN+V(I,Q)*SNV(I,Q)=-FM*SN+V(I,Q)*CNEND DODO I=1,NB(I)=A(I,I)END DOGOTO 10! 将特征值按大小排列15 M=N20 IF(M>0)THENJ=M-1M=0DO I=1,JIF(B(I).LT.B(I+1))THENB1=B(I)B(I)=B(I+1)B(I+1)=B1M=IDO K=1,NV1=V(K,I)V(K,I)=V(K,I+1)V(K,I+1)=V1END DOENDIFENDDOGOTO 20ENDIFEND5.3.6例对1991年5月1日至8月31日123天某纬圈上20个格点的850hPa高度场进行自然正交分解。
PROGRAM MAININTEGER,PARAMETER::P=20 !资料的空间点数INTEGER,PARAMETER::N=123 !资料的时间长度REAL(8),DIMENSION(P,N)::X,XF,ERRORREAL(8),DIMENSION(P)::XV !X的平均值INTEGER::LP=P;LW=3OPEN(12,FILE='F58Z850.DAT')DO IT=1,123READ(12,'(<LP>F8.1)')(X(I,IT),I=1,P)END DOCLOSE(12)XV=0DO I=1,PDO J=1,NXV(I)=XV(I)+X(I,J)END DOXV(I)=XV(I)/NEND DODO I=1,PDO J=1,NX(I,J)=X(I,J)-XV(I)END DOEND DOOPEN(12,FILE='EOF.DAT')CALL EOF(X,P,N,LW,XF)ERROR=X-XFDO I=1,PDO J=1,NXF(I,J)=XF(I,J)+XV(I)END DOEND DOWRITE(12,'(" 恢复场")')WRITE(12,'(<LP>F8.1)')((XF(I,J),J=1,P),I=1,N)WRITE(12,'(" 误差场")')WRITE(12,'(<LP>F8.1)')((ERROR(I,J),J=1,P),I=1,N)CLOSE(12)END输出结果为：特征值：1 2 3 4 5 6 7 8.2488D+06 .3194D+05 .1649D+05 .8222D+04 .4707D+04 .3137D+04 .1538D+04 .9665D+039 10 11 12 13 14 15 16.5258D+03 .4200D+03 .3926D+03 .3108D+03 .2749D+03 .2655D+03 .2395D+03 .1963D+0317 18 19 20.1922D+03 .1148D+03 .1110D+03 .9482D+02解释的方差(%)1 2 3 4 5 6 7 8 9 1078.01 88.02 93.19 95.77 97.25 98.23 98.71 99.02 99.18 99.3111 12 13 14 15 16 17 18 19 2099.44 99.53 99.62 99.70 99.78 99.84 99.90 99.94 99.97 100.00可见第一主分量解释的方差达78.01%，前三个主分量解释的方差达93.19%。