Math: 敏感度系数

Introduction

敏感度系数，又称灵敏度，表示项目评价指标对不确定因素的敏感程度。

计算公式：

\[E = \frac{{\Delta A}}{{\Delta F}}\]

式中，E为评价指标A对于不确定因素F的敏感度系数；

ΔA为不确定因素F发生ΔF变化率时，评价指标A的相应变化率，%；

ΔF为不确定因素F的变化率，%。

对于博士论文

情景情节设置变化率（r）使用在斜率因子（a），ΔF计算如下：

\[\begin{array}{c} \Delta F = \frac{{{I_{{\rm{changed}}}} - {I_{{\rm{no\_changed}}}}}}{{{I_{{\rm{no\_changed}}}}}}\\ = \frac{{a\left( {1 + r} \right)x + b - \left( {ax + b} \right)}}{{ax + b}}\\ = \frac{{ax + arx + b - ax - b}}{{ax + b}}\\ = \frac{{arx}}{{ax + b}} \end{array}\]

References

[1] 娄伟. 情景分析理论与方法[M]. 北京: 社会科学文献出版社, 2012, pp: 246~247.

[2] 敏感性分析详细过程.

Math: Moving Average

Introduction

In statistics, a moving average (rolling average or running average) is a calculation to analyze data points by creating series of averages of different subsets of the full data set. It is also called a moving mean (MM) or rolling mean and is a type of finite impulse response filter. Variations include: simple, and cumulative, or weighted forms (described below).

Given a series of numbers and a fixed subset size, the first element of the moving average is obtained by taking the average of the initial fixed subset of the number series. Then the subset is modified by "shifting forward"; that is, excluding the first number of the series and including the next number following the original subset in the series. This creates a new subset of numbers, which is averaged. This process is repeated over the entire data series. The plot line connecting all the (fixed) averages is the moving average. A moving average is a set of numbers, each of which is the average of the corresponding subset of a larger set of datum points. A moving average may also use unequal weights for each datum value in the subset to emphasize particular values in the subset.

已知一组序列和一个确定的子集大小，移动平均的第一个元素来自该序列的初始子集。而后子集继续向前移动，也就是，剔除序列的第一个成员再加入原序列的后一个成员。这就刚刚创建一个新的子集，求平均值。此过程在全序列不断重复至结束。图上与全部平均值关联的线段就是移动平均线。一个移动平均就是一个数值集合，每一元素对应较长数值序列的相应子集的平均值。移动平均可能采用加权计算，来突出特定样点在子集中的重要程度。

A moving average is commonly used with time series data to smooth out short-term fluctuations and highlight longer-term trends or cycles. The threshold between short-term and long-term depends on the application, and the parameters of the moving average will be set accordingly. For example, it is often used in technical analysis of financial data, like stock prices, returns or trading volumes. It is also used in economics to examine gross domestic product, employment or other macroeconomic time series. Mathematically, a moving average is a type of convolution and so it can be viewed as an example of a low-pass filter used in signal processing. When used with non-time series data, a moving average filters higher frequency components without any specific connection to time, although typically some kind of ordering is implied. Viewed simplistically it can be regarded as smoothing the data.

移动平均常应用在时间序列数据，平滑短期的波动并强调长期的趋势（周期）。短期和长期之间的阈值取决于实际情况，此时的移动平均参数须有对应的配置。例如，它常用在金融数据的技术分析，股价……数学上，一个移动平均就是一种卷积类型，所以他可以被视为信号处理的低通滤波器。当应用在非时间序列上，移动平均滤波不附带时间属性的高频率组分，即是该序列顺序有一定潜在意义。简单讲，移动平均是光滑数据的一种方式。

In financial applications a simple moving average (SMA) is the unweighted mean of the previous n data. However, in science and engineering the mean is normally taken from an equal number of data on either side of a central value. This ensures that variations in the mean are aligned with the variations in the data rather than being shifted in time. An example of a simple equally weighted running mean for a n-day sample of closing price is the mean of the previous n days' closing prices. If those prices are p_M, p_M-1, ..., p_M-(n-1) then the formula is

\[SMA = \frac{{{p_M} + {p_{M - 1}} + \cdots + {p_{M - \left( {n - 1} \right)}}}}{n}\]

还有Cumulative moving average，Weighted moving average，Exponential moving average，

时间尺度n确定方法

时间尺度n在不同的应用领域数值应是不同的，可以通过一定的算法估算具体情况下n的大小。移动平均（如简单移动平均等）并不一定有便捷的途径计算n的取值，但中心移动平均可以快速的求算n。所以，不妨利用中心移动平均求算n的途径为其他移动平均（如简单移动平均等）计算n的取值。

第一步，设置多种时间尺度（n=2,3……），原时间序列应用中心移动平均（Centered Moving Average）求得新的时间序列。

第二步，计算不同时间尺度下原时间序列与新输出时间序列的绝对差值代数和（Sum of absolute differences）。

第三步，图形展示不同时间尺度的代数和。曲线上的拐点（Knee Point of the Curve）即是对应于原时间序列的最佳时间尺度n。拐点的二阶导数等于0。

Fig. 1

如Fig. 1所示，(a)是原时间序列，(b)指出n等于9时是绝对差值代数和曲线的拐点，(c)黑线是原时间序列设定n为9时的中心移动平均曲线。

File.rar

References

[1] Moving average, from Wikipedia.

[2] Moving Average, Investopedia.

[3] How to decide window size for a moving average filter?

[4] Knee Point, from MathWorks File Exchange.

Matlab: A Practical Guide to Wavelet Analysis

Introduction

围绕C. Torrence and G. Compo (1998)文章示例辨析小波分析过程。

Data

These include the Nino3 sea surface temperature (SST) used as a measure of the amplitude of the El Nino-Southern Oscillation (ENSO). The Nino3 SST index is defined as the seasonal SST averaged over the central Pacific (5°S-5°N, 90°-150°W, Fig. 1). Data for 1871~1996 are from an area average of the U.K. Meteorological Office GISST2.3.

Fig. 1

示例数据如Fig. 2所示。

Fig. 2

Code Interpreter

归一化

1 2	variance = std(sst)^2; sst = (sst - mean(sst))/sqrt(variance) ;

首先，代码将原始数据做归一化处理，归一化并不是必要步骤，这里为方便比较而归一化，公式如下：

\[y = \frac{{x - m}}{s}\]

式中：y是归一化结果；m是x的平均值；s是x的标准差。

小波变换（Wavelet Transform，WT）

1 2 3 4 5 6 7 8 9 10

n = length(sst); % number of anomaly dt = 0.25 ; % sampling rate time = [0:length(sst)-1]*dt + 1871.0 ; % construct time array xlim = [1870,2000]; % plotting range pad = 1; % pad the time series with zeroes (recommended) dj = 0.25; % this will do 4 sub-octaves per octave s0 = 2*dt; % this says start at a scale of 6 months j1 = 7/dj; % this says do 7 powers-of-two with dj sub-octaves each lag1 = 0.72; % lag-1 autocorrelation for red noise background mother = 'Morlet';

参数设定。

File.rar

References

[1] Torrence, C.; Compo, G.P. A practical guide to wavelet analysis. Bull. Amer. Meteorol. Soc. 1998, 79, 61-78.

[2] Wavelet Software for Matlab.

Math: 时序多指标综合评价方法

Materials.

Math: 全局经济效益的综合评价方法

Materials.

Math: 基于标准差和平均差的权系数确定方法

File.rar.

Math: 客观赋权方法列表

熵权法：

\[{w_j} = \frac{{1 - k\sum\limits_{i = 1}^m {{f_{ij}}\ln {f_{ij}}} }}{{n - \sum\limits_{j = 1}^n {\left( {k\sum\limits_{i = 1}^m {{f_{ij}}\ln {f_{ij}}} } \right)} }},{f_{ij}} = \frac{{{r_{ij}}}}{{\sum\limits_{i = 1}^m {{r_{ij}}} }},k = \frac{1}{{\ln m}}\]

变异系数法：

\[{w_j} = \frac{{{\raise0.7ex\hbox{${\sqrt {\frac{1}{m}\sum\limits_{i = 1}^m {{{\left( {{r_{ij}} - \overline {{r_{ij}}} } \right)}^{\left( 2 \right)}}} } }$} \!\mathord{\left/ {\vphantom {{\sqrt {\frac{1}{m}\sum\limits_{i = 1}^m {{{\left( {{r_{ij}} - \overline {{r_{ij}}} } \right)}^{\left( 2 \right)}}} } } {\left| {\overline {{r_j}} } \right|}}}\right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{${\left| {\overline {{r_j}} } \right|}$}}}}{{\sum\limits_{j = 1}^n {\left( {{\raise0.7ex\hbox{${\sum\limits_{i = 1}^m {{{\left( {{r_{ij}} - \overline {{r_j}} } \right)}^{\left( 2 \right)}}} }$} \!\mathord{\left/ {\vphantom {{\sum\limits_{i = 1}^m {{{\left( {{r_{ij}} - \overline {{r_j}} } \right)}^{\left( 2 \right)}}} } {\left| {\overline {{r_j}} } \right|}}}\right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{${\left| {\overline {{r_j}} } \right|}$}}} \right)} }},\overline {{r_j}} = \frac{1}{m}\sum\limits_{i = 1}^m {{r_{ij}}} \]

均方差法

\[{w_j} = \frac{{\sqrt {\frac{1}{m}\sum\limits_{i = 1}^m {{{\left( {{r_{ij}} - \overline {{r_j}} } \right)}^2}} } }}{{\sum\limits_{j = 1}^n {\sqrt {\frac{1}{m}\sum\limits_{i = 1}^m {{{\left( {{r_{ij}} - \overline {{r_j}} } \right)}^2}} } } }},\overline {{r_j}} = \frac{1}{m}\sum\limits_{i = 1}^m {{r_{ij}}} \]

离差最大化法

\[{w_j} = {\textstyle{{\sum\limits_{i = 1}^m {\sum\limits_{k = 1}^m {\left| {{r_{ij}} - {r_{kj}}} \right|} } } \over {\sum\limits_{j = 1}^m {\sum\limits_{i = 1}^n {\sum\limits_{k = 1}^n {\left| {{r_{ij}} - {r_{kj}}} \right|} } } }}}\]

复相关系数法

\[{w_j} = \frac{{{\raise0.7ex\hbox{$1$} \!\mathord{\left/ {\vphantom {1 {{p_j}}}}\right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{${{p_j}}$}}}}{{\sum\limits_{j = 1}^n {\left( {{\raise0.7ex\hbox{$1$} \!\mathord{\left/ {\vphantom {1 {{p_j}}}}\right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{${{p_j}}$}}} \right)} }}\]

式中：p_j为x_j对x₁……x_n的复相关系数。

References

[1] 倪广亚, 刘学录, 李沁汶, 等. 基于数据信息特征的土地资源评价客观赋权方法的研究. 中国农学通报, 2014, 30(20): 255~262.

Math: The Bivariate Normal Distribution

Introduction

设（X, Y）为二维随机变量，若D（X），D（Y），Cov（X, Y）存在，且D（X），D（Y）都大于0，则有

\[{\rho _{XY}} = \frac{{Cov\left( {X,{\rm{ }}Y} \right)}}{{\sqrt {D\left( X \right)} \cdot \sqrt {D\left( Y \right)} }}\]

ρ_XY为相关系数。

假设随机变量X、Y的相关系数ρ_XY存在。当二者相互独立时，则E（XY）=E（X）E（Y），此时Cov（X, Y）=E（XY）—E（X）E（Y）=0，从而ρ_XY=0，即X、Y不相关。反之，若X、Y不相关，X、Y却不一定相互独立。其实，从“不相关”和“相互独立”的含义来看是明显的，不相关只是就线性关系来说的，而相互独立是就一般关系而言的。X、Y相互独立意味着两个变量之间没有任何关系，而X、Y不相关，仅仅说明X、Y之间无线性关系，但并不排除有非线性关系，比如对数关系、平方关系等等。因此，“不相关”是一个比“相互独立”弱得多的概念。

相关系数

设（X, Y）~N（μ₁, μ₂, σ₁, σ₂, ρ），它的概率密度为

\[\begin{array}{c} f\left( {x,y} \right) = \frac{1}{{2\pi {\sigma _1}{\sigma _2}\sqrt {1 - {\rho ^2}} }}\\ \cdot \exp \left\{ {\frac{{ - 1}}{{2\left( {1 - {\rho ^2}} \right)}}\left[ {\frac{{{{\left( {x - {\mu _1}} \right)}^2}}}{{\sigma _1^2}} - 2\rho \frac{{\left( {x - {\mu _1}} \right)\left( {y - {\mu _2}} \right)}}{{{\sigma _1}{\sigma _2}}} + \frac{{{{\left( {y - {\mu _2}} \right)}^2}}}{{\sigma _2^2}}} \right]} \right\} \end{array}\]

试求X和Y的相关系数。

解：可知

\[E\left( X \right) = {\mu _1},E\left( Y \right) = {\mu _2},D\left( X \right) = \sigma _1^2,D\left( Y \right) = \sigma _2^2\]

而

\[\begin{array}{c} {\rm{Cov}}\left( {X,Y} \right) = \int_{ - \infty }^{ + \infty } {\int_{ - \infty }^{ + \infty } {\left( {x - {\mu _1}} \right)\left( {y - {\mu _2}} \right)f\left( {x,y} \right){\rm{d}}} } x{\rm{d}}y\\ \left( {{\rm{Let }}t = \frac{1}{{\sqrt {1 - {\rho ^2}} }}\left( {\frac{{y - {\mu _2}}}{{{\sigma _2}}} - \rho \cdot \frac{{x - {\mu _1}}}{{{\sigma _1}}}} \right),u = \frac{{x - {\mu _1}}}{{{\sigma _1}}}} \right)\\ = \frac{1}{{2\pi }}\int_{ - \infty }^{ + \infty } {\int_{ - \infty }^{ + \infty } {\left( {{\sigma _1}{\sigma _2}\sqrt {1 - {\rho ^2}} tu + \rho {\sigma _1}{\sigma _2}{u^2}} \right){e^{ - \frac{{{u^2} + {t^2}}}{2}}}{\rm{d}}t{\rm{d}}u} } \\ = \frac{{\rho {\sigma _1}{\sigma _2}}}{{2\pi }}\left( {\int_{ - \infty }^{ + \infty } {{u^2}{e^{ - \frac{{{u^2}}}{2}}}{\rm{d}}u} } \right)\left( {\int_{ - \infty }^{ + \infty } {{e^{ - \frac{{{t^2}}}{2}}}{\rm{d}}t} } \right) + \frac{{{\sigma _1}{\sigma _2}\sqrt {1 - {\rho ^2}} }}{{2\pi }}\left( {\int_{ - \infty }^{ + \infty } {u{e^{ - \frac{{{u^2}}}{2}}}{\rm{d}}u} } \right)\left( {\int_{ - \infty }^{ + \infty } {t{e^{ - \frac{{{t^2}}}{2}}}{\rm{d}}t} } \right)\\ = \frac{{\rho {\sigma _1}{\sigma _2}}}{{2\pi }} \cdot \sqrt {2\pi } \cdot \sqrt {2\pi } \end{array}\]

即有

\[{\rm{Cov}}\left( {X,Y} \right) = \rho {\sigma _1}{\sigma _2}\]

于是

\[{\rho _{XY}} = \frac{{{\rm{Cov}}\left( {X,Y} \right)}}{{\sqrt {D\left( X \right)} \sqrt {D\left( Y \right)} }} = \rho \]

这就是说，二维正态随机变量（X，Y）的概率密度中的参数ρ就是X和Y的相关系数，因而二维正态随机变量的分布完全可由X、Y各自的期望、方差以及他们的相关系数确定。

References

[1] 李正耀, 周德强. 大学数学——概率论与数理统计. 北京: 科学出版社. 2009. Pp: 98-101.

Matlab: 小波分析—时间序列的多时间尺度分析（二）

一则小波分析帖子在我的科学网博客上获得了非常大的关注，至今这一帖子保持着留言询问最多的地位。不过，由于2015年11月的一场浩劫，这个帖子附带的示例代码就都烟消云散了，其后，还是有很多同学留言询问，我在这里把代码再回顾一遍，文字基本不变。

本部分对应——Matlab: 小波分析—时间序列的多时间尺度分析（2）。

Introduction

从原帖第3步计算小波系数开始重新讨论，与原帖不同，这里将尝试调用函数cwt()将最大尺度上限延伸至覆盖整个时间序列，得到小波系数，而后依次绘制、分析多时间尺度下的小波波动情况。

Fig. 1

Fig. 1为小波系数实部等值线图。这里呈现约为40年的主振荡周期，在整个时间序列上出现2个偏多中心和1个偏少中心，分别在1909、1998年和1950年。

Fig. 2

小波系数的模值是不同时间尺度变化周期所对应的能量密度在时间域中分布的反映，系数模值愈大，表明其所对应时段或尺度的周期性就愈强。从图 2可以看出，在降水量演化过程中，118~128年时间尺度模值最大（大于1000），横贯整个时间序列。

Fig. 3

小波系数的模方相当于小波能量谱，它可以分析出不同周期的震荡能量。由Fig. 3可知，118~128年时间尺度的能量最强、周期最显著，占据整个研究时域（1894~2010年）。

Fig. 4

小波方差图中（Fig. 4）存在较为明显的峰值。其中，最大峰值对应着78a的时间尺度，说明78a左右的周期震荡最强，为年降水量变化的第一主周期；55a时间尺度对应着第二峰值，为第二主周期（与原帖相同）；第三、第四峰值分别对应着30a和19a的时间尺度，它们依次为降水量的第三和第四主周期。这说明上述4个周期波动控制着降水量在整个时间域内的变化特征。

值得一提的是，小波方差在100~128a时间尺度上一直处于上升趋势，未出现拐点，这一现象值得再做考虑。

根据小波方差检验的结果，绘制演变的第一和第二主周期小波系数图。

Fig. 5

Fig. 6

从主周期趋势图中可以分析出在不同的时间尺度下，降水量存在的平均周期及丰—枯变化特征。Fig. 5显示，在78a特征时间尺度上，降水量变化的平均周期为50a左右，大约经历了2个丰—枯转换期；而在55a特征时间尺度上（Fig. 6），平均变化周期为35a左右，大约3个周期的丰—枯变化。

附上本帖数据及代码.

More

[1] A Practical Guide to Wavelet Analysis.

[2] 初学应用举例（小波分析，实部画图，小波方差画图及小波模画图等）.

[3] Matlab: 小波分析—时间序列的多时间尺度分析（一）.

Matlab: 小波分析—时间序列的多时间尺度分析（一）

一则小波分析帖子在我的科学网博客上获得了非常大的关注，至今这一帖子保持着留言询问最多的地位。不过，由于2015年11月的一场浩劫，这个帖子附带的示例代码就都烟消云散了，其后，还是有很多同学留言询问，我在这里把代码再回顾一遍，文字基本不变。

本部分对应——MATLAB：小波分析—时间序列的多时间尺度分析。

Introduction

时间序列（Time Series）是地学研究中经常遇到的问题。在时间序列研究中，时域和频域是常用的两种基本形式。其中，时域分析具有时间定位能力，但无法得到关于时间序列变化的更多信息；频域分析（如Fourier变换）虽具有准确的频率定位功能，但仅适合平稳时间序列分析。然而，地学中许多现象（如河川径流、地震波、暴雨、洪水等）随时间的变化往往受到多种因素的综合影响，大都属于非平稳序列，它们不但具有趋势性、周期性等特征，还存在随机性、突变性以及“多时间尺度”结构，具有多层次演变规律。对于这类非平稳时间序列的研究，通常需要某一频段对应的时间信息，或某一时段的频域信息。显然，时域分析和频域分析对此均无能为力。

20世纪80年代初，由Morlet提出的一种具有时—频多分辨功能的小波分析（Wavelet Analysis）为更好的研究时间序列问题提供了可能，它能清晰的揭示出隐藏在时间序列中的多种变化周期，充分反映系统在不同时间尺度中的变化趋势，并能对系统未来发展趋势进行定性估计。

本帖以美国某气象站1894~2010年连续的年降水量为例，试用小波分析，完成如下任务：①小波变换系数；②绘制小波系数实部等值线图；③绘制小波系数模和模方等值线图；④绘制小波方差图；以及⑤绘制不同时间尺度的小波实部过程线。所谓年降水量时间序列的多时间尺度是指：年降水量在演化过程中，并不存在真正意义上的变化周期，而是其变化周期随着研究尺度的不同而发生相应的变化，这种变化一般表现为小时间尺度的变化周期往往嵌套在大尺度的变化周期之中。也就是说，年降水量变化在时间域中存在多层次的时间尺度结构和局部变化特征。

小波分析的计算过程请参考：小波分析—经典、小波方差制作步骤等。

小波分析的基本理论在此不多叙述，请参考其他文献。本帖主要介绍小波分析的一般过程，数据及部分代码将附在文末。

0．年降水量的变化趋势分析

该站点的年降水量变化情况（LI_plot函数）如图 0所示，其发展呈微微上升的趋势，降水量最高年份出现在2003年，全年累计达到1610.70mm，最低值出现在1965年，累计仅有748.90mm。

Fig. 0

1．数据的加载

加载多年降水量至MATLAB使用Load命令。

2．边界效应的消除或减小

由于本例中的实测年降水量数据为有限时间数据序列，在时间序列的两端可能会产生“边界效用”。为消除或减小序列开始点和结束点附近的边界效应，须对其两端数据进行延伸。在进行完小波变换后，去掉两端延伸数据的小波变换系数，保留原数据序列时段内的小波系数。本例中，我们利用Matlab小波工具箱中的信号延伸（Signal Extension）功能，对数据两端进行对称性延伸。

具体步骤：在Matlab界面的“Command Window”中输入小波工具箱调用命令“wavemenu”，按Enter键弹“Wavelet Toolbox Main Menu”（小波工具箱主菜单）界面（Fig. 1上）；然后单击“Signal Extension”，打开Signal Extension /Truncation窗口，单击“File”菜单下的“Load Signal”，选择prec.mat文件单击“打开”，出现信号延伸界面。Matlab的Extension Mode菜单下包含多种延伸方式和Direction to extend菜单下的3种延伸模式（Both、Left and Right），在这里我们选择对称性两端延伸进行计算。具体操作过程是：在Extension Mode下选择“Symmetric(Whole-Point)”，Dircetion to extend下选择“Both”，单击“Extend”按钮进行对称性两端延伸计算（Fig. 1下），然后单击“File”菜单下的“Save Tranformed Signal”，将延伸后的数据结果存为ePrec.mat文件。

Fig. 1

从ePrec文件可知（entendvalue函数），系统自动将原时间序列数据向前对称延伸5个单位，向后延伸6个单位。

3．计算小波系数

选择Matlab小波工具箱中的Morlet复小波函数对延伸后的数据序列（ePrec.mat）进行小波变换，计算小波系数并保存。

小波工具箱主菜单界面见Fig. 1上，单击“One-Dimensional”下的子菜单“Complex Continuous Wavelet 1-D”，打开一维复连续小波界面，单击“File”菜单下的“Load Signal”按钮，载入时间序列数据ePrec.mat。Fig. 2第一排的左侧为信号显示区域，右侧区域给出了信号序列和复小波变换的有关信息和参数，主要包括数据长度（Data Size）、小波函数类型（Wavelet：cgau、shan、fbsp和cmor）、取样周期（Sampling Period）、周期设置（Scale Setting）和运行按钮（Analyze），以及显示区域的相关显示设置按钮。本例中，我们选择cmor (1-1.5)、取样周期为1、最大尺度为64（延伸后的时间序列的一半），单击“Analyze”运行按钮，计算小波系数。然后单击“File”菜单下的“Save Coefficients”，保存小波系数为cePrec.mat文件。

Fig. 2

4．计算Morlet复小波系数的实部

去除两端延伸数据的小波系数（entendvalue函数），并计算小波系数实部（real函数）。

5．绘制小波系数实部等值线图

这部分过程应用LI_contourf函数，年降水量小波系数实部等值线图见Fig. 3。

如Fig. 3所示的小波系数实部等值线图。其中，横坐标为时间（年份），纵坐标为时间尺度，图中的等值曲线为小波系数实部值。小波系数实部等值线图能反映年降水量序列不同时间尺度的周期变化及其在时间域中的分布，进而能判断在不同时间尺度上，年降水量的未来变化趋势。

Fig. 3

由Fig. 3可以看出降水量1894~2010年演化过程中存在着28~36a的主振荡周期（这一周期是多个振荡周期叠加的矢量和，随后将在小波方差图中对这一主周期进行分解，剥离出第一主周期等等）。在整个时间尺度上出现4个偏多中心和3个偏少中心，分别为1899、1938、1972、2007年和1919、1956、1990年。

红框处出现4个偏多中心和3个偏少中心，一对临近的偏多和偏少波动就可以被识别为一个周期，所以是28~36a主震荡周期。

周期判断可以从正弦曲线上类比，如下图，两个最邻近波峰（或波谷）之间是一个完整的周期，或者两个最邻近的与x轴向下（或向上）交点也是一个完整的周期。

那么按照正弦曲线的周期检测思路，两个最邻近的偏多中心（或偏少中心）就是一个完整的周期长度。

File.rar.

6．绘制小波系数模和模方等值线图

首先，计算小波系数的模和模方，再利用LI_contourf函数绘制模和模方等值线图。

Fig. 4

Morlet小波系数的模值是不同时间尺度变化周期所对应的能量密度在时间域中分布的反映，系数模值愈大，表明其所对应时段或尺度的周期性就愈强。从Fig. 4可以看出，在降水量演化过程中，纵轴40~64a时间尺度模值最大（大于400），但1924~1994年之间模值小余200，说明在此时段内40~64年时间尺度的周期变化并不明显，1995年之后模值再次增大，此时期后40~64年时间尺度的周期变化趋于显著。

Fig. 5

小波系数的模方相当于小波能量谱，它可以分析出不同周期的震荡能量。由Fig. 5可知，40~64年时间尺度的能量最强、周期最显著，但它的周期变化具有局部性（1924年之前和1994年之后）；10~15年时间尺度能量虽然较弱，但周期分布比较明显，几乎占据整个研究时域（1894~2010年）。

7．绘制小波方差图

小波方差的计算过程参考：小波方差制作步骤（参考文献存在一个错误，小波方差未除N，本贴小波方差计算已有修正）。

小波方差图能反映降水量时间序列的波动能量随尺度（a）的分布情况。可用来确定降水量演化过程中存在的主周期。

Fig. 6

8．主周期趋势图的绘制及其在多时间尺度分析中的作用

根据小波方差检验的结果，绘制演变的第一和第二主周期小波系数图。

Fig. 7

Fig. 8

从主周期趋势图中可以分析出在不同的时间尺度下，降水量存在的平均周期及丰-枯变化特征。Fig. 7显示，在55a特征时间尺度上，降水量变化的平均周期为35a左右，大约经历了3个丰—枯转换期；而在30a特征时间尺度上（Fig. 8），平均变化周期为20a左右，大约6个周期的丰—枯变化。

File.rar.

More

[1] A Practical Guide to Wavelet Analysis.

[2] 初学应用举例（小波分析，实部画图，小波方差画图及小波模画图等）.

Math: DEA

File.rar.

Math: 可拓综合评价法的建模步骤

矛盾问题，就是指在现有条件下无法实现人们要达到的目标的问题。

Summary

利用可拓综合评价法，一般要经历确定经典域、确定节域、确定待评价物元、确定评价指标的权重、确定待评价事物关于各类别的关联度以及确定待评价事物的类别和级别变量特征值六个步骤。

Step 1 确定经典域

设有m个评价物元（或评价类别）N₁，N₂，…，N_m，将各评价物元（或类别）对应的特征值范围用[a_ij, b_ij]表示，则同征物元R₀可表示为：

\[{R_0} = \left[ {\begin{array}{*{20}{c}} N&{{N_1}}&{{N_2}}& \cdots &{{N_m}}\\ {{c_1}}&{\left[ {{a_{11}},{b_{11}}} \right]}&{\left[ {{a_{12}},{b_{12}}} \right]}& \cdots &{\left[ {{a_{1m}},{b_{1m}}} \right]}\\ {{c_2}}&{\left[ {{a_{21}},{b_{21}}} \right]}&{\left[ {{a_{22}},{b_{22}}} \right]}& \cdots &{\left[ {{a_{2m}},{b_{2m}}} \right]}\\ \vdots & \vdots & \vdots & \vdots & \vdots \\ {{c_n}}&{\left[ {{a_{n1}},{b_{n1}}} \right]}&{\left[ {{a_{n2}},{b_{n2}}} \right]}& \cdots &{\left[ {{a_{nm}},{b_{nm}}} \right]} \end{array}} \right]\]

式中：N_j——第j个评价物元（或类别）；c_i——第i个评价指标；v_ij=[a_ij,b_ij]——N_j关于c_i所规定的量值范围，即经典域。

Step 2 确定节域

\[{R_p} = \left( {P,C,{V_p}} \right) = \left[ {\begin{array}{*{20}{c}} P&{{c_1}}&{\left[ {{a_{1p}},{b_{1p}}} \right]}\\ {}&{{c_2}}&{\left[ {{a_{2p}},{b_{2p}}} \right]}\\ {}& \vdots & \vdots \\ {}&{{c_n}}&{\left[ {{a_{np}},{b_{np}}} \right]} \end{array}} \right]\]

式中：P——评价类别的全体；[a_ip，b_ip]——P关于c_i所取的量值范围，即节域。

Step 3 确定待评价物元

对待评价物元，把监测数据或分析结果用物元表示为：

\[{R_d} = \left[ {\begin{array}{*{20}{c}} P&{{c_1}}&{{v_1}}\\ {}&{{c_2}}&{{v_2}}\\ {}& \vdots & \vdots \\ {}&{{c_n}}&{{v_n}} \end{array}} \right]\]

式中：R_d——待评价物元；v_i——待评价事物对应于c_i的数值。

Step 4 确定评价指标的权重

\[{W_i} \ge 0\left( {i = 1,2, \cdots ,n} \right)\] \[\sum\limits_{i = 1}^n {{W_i}} = 1\]

Step 5 确定关联度

计算距：

\[\begin{array}{c} \rho \left( {{v_i},{v_{ij}}} \right) = \left| {{v_i} - \frac{1}{2}\left( {{a_{ij}} + {b_{ij}}} \right)} \right| - \frac{1}{2}\left( {{b_{ij}} - {a_{ij}}} \right)\\ \rho \left( {{v_i},{v_{ip}}} \right) = \left| {{v_i} - \frac{1}{2}\left( {{a_{ip}} + {b_{ip}}} \right)} \right| - \frac{1}{2}\left( {{b_{ip}} - {a_{ip}}} \right) \end{array}\]

式中：ρ（v_i,v_ij）——点v_i与区间v_ij的距；ρ（v_i,v_ip）——点v_i与区间v_ip的距。

计算关联函数：

\[{K_j}\left( {{v_i}} \right) = \left\{ {\begin{array}{*{20}{c}} {\frac{{\rho \left( {{v_i},{v_{ij}}} \right)}}{{\rho \left( {{v_i},{v_{ip}}} \right) - \rho \left( {{v_i},{v_{ij}}} \right)}}{v_i} \notin {v_{ij}}}\\ {\frac{{ - \rho \left( {{v_i},{v_{ij}}} \right)}}{{\left| {{v_{ij}}} \right|}}{v_i} \in {v_{ij}}} \end{array}} \right.\]

K_j（v_i）——关联函数，待评价事物的指标c_i关于类别j的归属度；|v_ij|——区间[a_ij，b_ji]的长度，即|b_ij-a_ij|。

计算关联度：

\[{K_j}\left( p \right) = \sum\limits_{i = 1}^n {{W_i}{K_j}\left( {{v_i}} \right)} \]

式中：K_j（p）——在考虑指标权重下，待评价事物各指标c_i关于类别j的关联度的组合值。

Step 6 确定待评价事物的类别和级别变量特征值

若K_j0（p）=maxK_j（p）（j=1，2，…，m），则评定p属于类别j0。记：

\[\overline {{K_j}} \left( p \right) = \frac{{{K_j}\left( p \right) - \mathop {\min }\limits_{1 \le j \le m} {K_j}\left( p \right)}}{{\mathop {\max }\limits_{1 \le j \le m} {K_j}\left( p \right) - \mathop {\min }\limits_{1 \le j \le m} {K_j}\left( p \right)}}\]

则p的级别变量特征值j*为：

\[{j^ * } = \frac{{\sum\limits_{j = 1}^m {j \cdot \overline {{K_j}} \left( p \right)} }}{{\sum\limits_{j = 1}^m {\overline {{K_j}} \left( p \right)} }}\]

File.rar.

References

[1] 何逢标. 综合评价方法MATLAB实现. 北京: 中国社会科学出版社. pp: 269~270.

[2] 方统中, 杜耘, 蔡述明, 等. 模糊数学在洪湖营养化评价中的应用. 浙江林学院学报, 2008, 25(4):517~521.

[3] 高军省, 吕小凡. 湖泊富营养化评价的可拓学方法及应用. 长江大学学报（自然科学版）. 2010, 7(1):33~37.

Math: 可拓综合评价法概述

Summary

模糊综合评价法是建立在模糊集合基础上的评价方法，利用它可解决“内涵明确，外延不明确”的具有模糊性的一类评价问题。而可拓综合评价法是建立在可拓集合基础上的评价方法，它不仅可以从数量上刻画被评价对象本身在状态的所属程度，而且还可以从数量上刻画何时为一种性态与另一种性态的边界。

在日常生活中，我们常说“量变引起质变”，但究竟多少“量”的累积才会引起“质”的变化，通常难以说明。“可拓学是一种把事物的质与量有机结合起来的理论，它从定性与定量两个角度去研究解决具有矛盾的问题，其理论支柱是物元理论和可拓集合论，其逻辑细胞是物元”。

可拓学是我国学者蔡文于1983年创立的一门新学科，它用形式化的模型研究事物拓展的可能性和开拓创新的规律与方法。物元以有序的三元组R=（N，C，V）来表达。若事物N有n个特征C（c₁，c₂，…，c_n），待评价的m个物元R₁=（N₁，C，V₁），R₂=（N₂，C，V₂），…，R_m=（N_m，C，V_m）具有相同的特征C，则具有相同特征的物元体R可表示为：

\[R = \left[ {\begin{array}{*{20}{c}} N&{{N_1}}&{{N_2}}& \cdots &{{N_m}}\\ {{c_1}}&{{v_{11}}}&{{v_{12}}}& \cdots &{{v_{1m}}}\\ {{c_2}}&{{v_{21}}}&{{v_{22}}}& \cdots &{{v_{2m}}}\\ \vdots & \vdots & \vdots & \vdots & \vdots \\ {{c_n}}&{{v_{n1}}}&{{v_{n2}}}& \cdots &{{v_{nm}}} \end{array}} \right]\]

其中，N_i——待评价的物元；N——待评价物元N₁，N₂，…，N_m的全体；v_ij——第j个待评价物元第i个特征的量值，i=1，2，…，n；j=1，2，…，m。

以此为基础，建立的可拓综合评价法，通过具有相同特征的物元的表示、各特征量值范围的确定、待评价物元各项特征值的描述、关联度的计算，最终得出评价物元所属的类别，现今在水利、环境、经济与管理等众多领域得到了越来越多的应用。

References

[1] 何逢标. 综合评价方法MATLAB实现. 北京: 中国社会科学出版社. pp: 269~270.

Math: 熵权法

Description

物理学上的熵（entropy）用来表示任何一种能量在空间中的分布的均匀程度，能量分布的越均匀，熵就越大。当一个体系的能量完全均匀分布时，这个系统的熵就达到最大值。在信息论中，熵又称为平均信息量，它是信息的一个度量。

利用熵的概念确定指标权重的方法称为熵权法。其出发点是根据某同一指标观测值之间的差异程度来反映其重要程度，如果各被评价对象的某项指标的数据差异不大，则反映该指标对评价系统所起的作用不大。

熵权法是一种客观的赋权方法，它是利用各指标的熵值所提供的信息量的大小来决定指标权重的方法。熵权法的作用有：用熵权法给指标赋权可以避免各评价指标权重的人为因素干扰，使评价结果更符合实际；通过对各指标熵的计算，可以衡量出指标信息量的大小，从而确保所建立的指标能反映绝大部分的原始信息。

Method

（1）形成决策矩阵

设参与评价的对象集为M=（M₁，M₂，…，M_m），指标集为D=（D₁，D₂，…，D_n），评价对象M_i对指标D_j的值记为x_ij（i=1，2，…，m；j=1，2，…，n），则形成的决策矩阵X为：

\[X = \left[ {\begin{array}{*{20}{c}} {}&{{D_1}}&{{D_2}}& \cdots &{{D_n}}\\ {{M_1}}&{{x_{11}}}&{{x_{12}}}& \cdots &{{x_{1n}}}\\ {{M_2}}&{{x_{21}}}&{x{}_{22}}& \cdots &{{x_{2n}}}\\ \vdots & \vdots & \vdots & \vdots & \vdots \\ {{M_m}}&{{x_{m1}}}&{{x_{m2}}}& \cdots &{{x_{mn}}} \end{array}} \right]\]

（2）标准化决策矩阵

为了消除各指标量纲不同对方案决策带来的影响，或者处理一些指标值为负的决策问题，对决策矩阵X进行标准化处理，从而形成标准化矩阵V=（v_ij）_m×n。根据指标的性质，将指标分为两类。一类是越大越优型指标，也称为效益型指标；另一类是越小越优型指标，也称为成本型指标。

标准化处理时根据指标性质，采用相应的标准化形式：

对于越大越优型指标：

\[{v_{ij}} = \frac{{{x_{ij}} - \min \left( {{x_j}} \right)}}{{\max \left( {{x_j}} \right) - \min \left( {{x_j}} \right)}}\]

对于越小越优型指标：

\[{v_{ij}} = \frac{{\max \left( {{x_j}} \right) - {x_{ij}}}}{{\max \left( {{x_j}} \right) - \min \left( {{x_j}} \right)}}\]

（3）计算第j项指标下，第i个评价对象的特征比重

对于某一个指标j，v_ij的值差异越大，表明该项指标对于被评价对象的作用越大，即该项指标提供给被评价对象的有用信息越多。根据熵的概念，信息的增加意味着熵的减少，熵可以用来度量这种信息量的大小。

记第j项指标下，第i个评价对象的特征比重为p_ij，则

\[{p_{ij}} = \frac{{{v_{ij}}}}{{\sum\limits_{i = 1}^m {{v_{ij}}} }}\]

（4）计算第j项指标的熵值e_j

\[{e_j} = - \frac{1}{{\ln \left( m \right)}}\sum\limits_{i = 1}^m {{p_{ij}}\ln \left( {{p_{ij}}} \right)} \] 当p_ij=0或者p_ij=1时，认为p_ijln(p_ij)=0。

（5）计算第j项指标的差异性系数d_j

观察熵值的计算公式，对于某一项指标D_j，v_ij的差异越小，e_j越大。当各被评价对象第j项指标值全部相等时，e_j=e_max=1。根据熵的概念，各被评价对象第j项指标值差异越大，表明该指标反映的信息量越大。因此，定义差异系数d_j：

\[{d_j} = 1 - {e_j}\]

（6）确定各指标的熵权

\[{w_j} = \frac{{{d_j}}}{{\sum\limits_{k = 1}^n {{d_k}} }}j = 1,2, \cdots ,n\]

Parctices.rar.

References

[1] 郭亚军. 综合评价理论、方法及应用. 北京: 科学出版社, 2007. Pp15~16.

[2] 何鑫, 朱宏泉, 高成凤. 基于熵权法与TOPSIS法的房地产项目投资风险评价. 商业研究, 2009, 03:105~108.

Math: 评价指标的一致化

Description

一般来说，指标x₁, x₂, …, x_m中，可能含有“极大型”指标、“极小型”指标、“居中型”指标和“区间型”指标。对于某些定量指标，如产值、利润等，我们自然期望它们的取值越大越好，这类指标我们称之为极大型指标；而对于诸如成本、能耗等一类指标，我们自然期望它们的取值越小越好，这类指标称之为极小型指标；诸如人的身高、体重等指标，我们既不期望它们的取值越大越好，也不期望期望它们的取值越小越好，而是期望它们的取值越居中越好，我们称这类指标为居中型指标；而区间型指标是期望其取值以落在某个区间内为最佳的指标。

若指标x₁, x₂, …, x_m中既有极大型指标、极小型指标，又有居中型指标或区间型指标，则必须在对各备选方案进行综合评价之前，将评价指标的类型作一致化处理。否则，就无法定性地判定综合评价函数中的yi值是否是取值越大越好，或是取值越小越好，或是取值越居中越好。因此，也就无法根据y值的大小来综合评价各备选方案的优劣。因此，需将指标作类型一致化的处理。

对于极小型指标x，令

\[{x^*} = M - x\]

或

\[{x^*} = \frac{1}{x}\left( {x > 0} \right)\]

式中，M为指标x的一个允许上界。

对于居中型指标x，令

\[{x^*} = \left\{ {\begin{array}{*{20}{c}} {2\left( {x - m} \right),m \le x \le \frac{{M + m}}{2}}\\ {2\left( {M - x} \right),\frac{{M + m}}{2} \le x \le M} \end{array}} \right.\]

式中，m为指标x的一个允许下界，M为指标x的一个允许上界。

对于区间型指标x，令

\[{x^*} = \left\{ {\begin{array}{*{20}{c}} {1.0 - \frac{{{q_1} - x}}{{\max \left\{ {{q_1} - m,M - {q_2}} \right\}}},x < {q_1}}\\ {1.0,x \in \left[ {{q_1},{q_2}} \right]}\\ {1.0 - \frac{{x - {q_2}}}{{\max \left\{ {{q_1} - m,M - {q_2}} \right\}}},x > {q_2}} \end{array}} \right.\]

式中，[q₁,q₂]为指标x的最佳稳定区间，M、m分别为x的允许上、下界。

如此，非极大型评价指标x通过上式可以转换为极大型指标了。

References

[1] 郭亚军. 综合评价理论、方法及应用. 北京: 科学出版社, 2007. Pp15~16.

Math: Symbols for Preference Relations

Summary

≻succeeds和≺precedes的意义。

Set Theory Symbols.

References

[1] ≻.

[2] ≺.

[3] What does the math symbol ≺ mean, as in 1≺loglogn≺logn?

[4] What's the meaning of ≻≻ operator?

Math: Tips

[1] 绝对变化与相对变化？