首页 / 技术文章 / 使用 MATLAB^® 工具箱

使用 NAG MATLAB^® 工具箱 - Part 3

本系列的文章可参考 http://cn.nag-gc.com/IndustryArticles/usingtoolboxmatlab.asp 以及 http://cn.nag-gc.com/IndustryArticles/usingtoolboxmatlabpart2.asp；本文中，我们将持续探索我们的工具箱， 并演示如何在 MATLAB 中使用我们的 NAG 算法库，并利用 MATLAB 的图形功能显示结果。

注意： 本文中的范例皆是由我们的演示程序中截取的，若您直接以文中的 MATLAB 程序运行，也许无法完整的获得如文中的结果。文中范例的完整的程序，请您自 此处下载。

G13 函数 - 时间序列分析

所谓时间序列指的是在某段时间中，收集潜在时间相依流程中的观测点集合。NAG 的 G13 函数系列包含了许多调查与建构时间序列的统计模型； 由这些函数建构出的模型可以让我们更了解资料的特性，或者是从时间序列中进行预测 (也就是对未来的行为进行预测)。 例如：一般说的 自我回归集成移动平均 (autoregressive integrated moving average) (ARIMA) 模型能够预测时间序列，请看以下例子。

了解时间序列特性的第一步，是去计算它的 自相关函数，它描述了在不同时间点所存在的相关性 (依赖程度)。 将两个不同时间段进行分割称之为 lag，而自相关函数通常是依不同的 lag 值，以 自相关系数 来表示结果。 g13ab 函数能够依许多基本的统计量，如平均值与变异数来计算这个系数。以下是程序片段：

% 时间序列数据 (来自太阳黑子观测值) x = [5; 11; 16; 23; 36; 58; 29; 20; 10; 8; 3; 0; 0; 2; 11; 27; 47; 63; 60; 39; 28; 26; 22; 11; 21; 40; 78; 122; 103; 73; 47; 35; 11; 5; 16; 34; 70; 81; 111; 101; 73; 40; 20; 16; 5; 11; 22; 40; 60; 80.9]; % nk 表示 lags 值，自相关函数所需 nk = int32(40); % 调用 NAG 函数 [xm, xv, coeff, stat, ifail] = g13ab(x, nk);

因为 lag 是一个不连续的变量，自相关比较适合用直方图表示 (或者称 autocorrelogram)，如下图：

图 1：计算时间序列的自相关函数

自相关函数包含了潜在时间相依流程的定量与定性讯息；本例中，当以 ARIMA 模型去拟合时间序列，可看出大约具有 11 周期的摆动。 autocorrelation 图形的形状可以显现出模型参数。图 1 中，可以了解到这个时间序列为 非周期性的，因此还需进一步探索。 如果相关系数在前面几个 lags 中的值很高，然后快速下降，便是我们一般称的 移动平均 (MA) 序列， 若结果以波形方式呈现，便为 自回归 (AR) 序列。大部分的例子中，通常需要使用 ARIMA 模型 (同时拥有 AR 与 MA) 去拟合序列的资料。

除了自相关函数，我们也可以自 偏自相关 (partial autocorrelation) 函数 的图形中，获得其他更多的讯息。 我们可以在以上的例子中以 g13ac 函数取代 g13ab 函数的调用。

本范例的 MATLAB 脚本位于 NAGToolboxDemos/Time_series_analysis/g13ab_demo.m，也可下载 压缩文件。

D01 函数 - 积分

在数值分析中，定积分中的一维或多维度的求值是非常普遍的。D01 函数提供了非常多种的算法求解； 个别函数的调用全然依积分的形式而定。如果函数的形式可以分析出来，如被积函数包含了奇异值 (singularities)，尤其是代数或对数型的样式，那么 d01aj 函数最为合适。 其他较为特殊的函数，例如：如果被积函数为震荡函数，可以采用 d01ak；已知点为不连续的被积函数，则使用 d01al 函数。

d01aj 函数采用的是 自适应 算法 - 意思是它会将要积分的函数，切分成许多的子区间， 并反复切割，直到满足所设定的精确度条件为止 (以下程序中以 epsabs 与 epsrel 变量来指定)。

以下程序调用 NAG 函数求解积分

epsabs = 0.0; epsrel = 0.000001; % lw 与 liw 是输出矩阵 w 与 iw 的矩阵大小 % lw 必须介于 800 与 2000 之间；liw 为 lw 的四分之一 lw = int32(1000); liw = lw/4; iw = zeros(liw, 'int32'); w = zeros(lw, 'double'); % 调用 NAG 函数并验证误差 [result, abserr, w, iw, ifail] = d01aj(func, a, b, epsabs, epsrel); if ifail ~= 0 disp([' Non-zero ifail after d01aj call ', num2str(ifail)]); return; end;

此程序中，d01aj 函数中所调用的 func 变量，指定的是 MATLAB m 档案的函数字符串名； 这个档案必须是一个函数，将会传回给定值的积分值。 以下是一个例子：

function [result] = myfunc(x) result = x^2*abs(sin(5*x));

(这是下图中所计算的被积函数)。除了近似的积分估计值外，d01aj 函数还会输出最后子区间的误差值。 图 2 显示出每一个积分子区间的结果以及其相对应的误差值。 从这个积分中也可以看出，在函数快速变化的区域，此算法会窄化子区间的分割；且与较大误差有关的子区间，其区间宽度也会较小 (函数接近于 0 的地方)。

图 2：积分求值估计

本范例的 MATLAB 脚本位于 NAGToolboxDemos/Quadrature/d01aj_demo.m， 也可下载 压缩文件。

G03 函数 - 多变量方法

多变量数据集乃是包含许多不同 目标 所观察到的 变量。在科学分析中有许多这样的范例数据，NAG 算法库中的 G03 函数可以利用这些数据集来进行分析。 举例说明，在英国的环境科学家，想要了解在他们所观测的，分布在英国与欧洲大陆间的 14 个区域中，进行观测所取样 300 个水鼠 (Arvicola 属) 头颅表现的 13 个特征中，分属于多少种的水鼠种类。 。自欧洲取得的资料已经归类为二大类 (terrestis 与 sapidus)；科学家的任务便要来决定英国所取得的资料是属于那一个种类的。

数据处理是以 14 个区域中所观测的 13 个变量的平均值计算开始。这可以视为 13 维度空间中的 14 个点， 一般我们会用 主成分分析 (principal components analysis) - g03aa 函数 来分析这样的资料。主成分分析是用来将原始数据维度简化的一个方法；只要这些数值可以足够表示原来的数据，那么这些数值就可以取代原来的数值。 然而对于这个田鼠的数据而言，考虑到前三个主成分只能解释原始资料 65% 的变异，显然是不适合的。

另一个可以应用在这个问题上的技术就是 metric scaling。第一步，便是建立一个 14 乘 14 的 相异矩阵 (dissimilarity matrix) ， 其中的每个元素就是其与原来 13 维度空间的距离。g03ea 函数可用来计算此相异矩阵。以下是程序片段：

% 每一个数据区域的名称 label = char( 'Surrey', 'Shropshire', 'Yorkshire', ... 'Perthshire', 'Aberdeen', 'Eilean Gamhna', 'Alps', ... 'Yugoslavia', 'Germany', 'Norway', 'Pyrenees I', ... 'Pyrenees II', 'North Spain', 'South Spain' ); % 田鼠的资料。14 列 (每一个区域)，13 行 (每一个变量) data0 = [ [48.5 89.2 7.9 42.1 92.1 100.0 100.0 ... 35.3 11.4 100.0 71.9 31.6 2.8]; [67.7 67.0 2.0 23.0 93.0 100.0 86.0 ... 44.0 14.0 99.0 97.0 31.0 17.0]; [51.7 84.8 0.0 31.3 88.2 100.0 94.1 ... 18.8 25.0 100.0 83.3 33.3 5.9]; [42.9 50.0 0.0 50.0 77.3 100.0 90.9 ... 36.4 59.1 100.0 100.0 38.9 0.0]; [18.1 79.6 4.1 44.9 79.6 100.0 77.6 ... 16.7 37.1 100.0 90.4 9.8 0.0]; [65.0 81.8 9.1 31.8 81.8 100.0 59.1 ... 20.0 30.0 100.0 100.0 5.0 9.1]; [57.1 76.2 21.4 38.1 66.7 97.6 14.3 ... 23.5 9.5 100.0 91.4 11.8 17.5]; [26.7 53.1 23.5 38.2 44.1 94.1 11.8 ... 11.8 18.2 100.0 94.9 12.5 5.9]; [38.5 67.9 17.9 21.4 82.1 100.0 60.7 ... 35.7 24.0 100.0 91.7 37.5 0.0]; [33.3 83.3 27.8 29.4 86.1 100.0 63.9 ... 53.8 18.8 100.0 83.3 8.3 34.3]; [47.6 92.9 26.7 10.0 36.7 100.0 50.0 ... 14.3 7.4 100.0 86.4 90.9 3.3]; [60.0 90.9 13.6 68.2 40.9 100.0 18.2 ... 100.0 5.0 80.0 90.0 50.0 0.0]; [53.8 88.1 7.1 33.3 88.1 100.0 19.0 ... 85.7 9.8 73.8 72.2 73.7 2.4]; [29.2 74.0 16.0 46.0 86.0 100.0 18.0 ... 88.0 16.3 72.0 80.4 69.6 4.0]; ]; % 使用标准转换，消除区间中不一致的变量差异的影响 data1 = asin(1.0 - 0.02*data0); update = 'I'; % 距离为欧几里得平方 dist = 'S'; % 对变数不做尺度转换 scal = 'U'; % 观察值数量 (收集数据的区域) nobs = int32(14); % 变量的数量 (每一区域所量测的数值) nvar = int32(13); % 哪一个变量并不会包含在距离计算中 isx = ones(1, nvar, 'int32'); % 每一个变量的缩放调整参数 (未使用，但是必须提供) s = ones(1, nvar, 'double'); % 输入距离矩阵 (未使用，必须提供) d = zeros(nobs*(nobs-1)/2, 1, 'double'); % 调用 NAG 函数计算相异矩阵 [sOut, dOut, ifail] = g03ea(update, dist, scal, nobs, data1, ... isx, s, d, 'm', nvar);

第二个步骤便是找出一个新的矩阵能够表示较少维度 (本例中是 3) 的数据距离；我们使用 metric scaling 将旧的矩阵投射到新的矩阵中。 接下来的 g03fa 函数，便是用来处理 metric scaling：

% 依变量数量将矩阵做标准化处理 dOut = dOut ./ double(nvar); % 每个区域中有不同观察数量，修正偏差 ns = [19 50 17 11 49 11 21 17 14 18 16 11 21 25]; jend = nobs-1; for j = 1:jend istart = j+1; for i = istart:nobs index = (i-1)*(i-2)/2 + j; dOut(index) = abs(dOut(index) - (1.0/ns(i) + 1.0/ns(j))); end end % 进行 nvar 到 ndim 的 metric scaling 计算 roots = 'L'; ndim = int32(3); [x, eval, ifail] = g03fa(roots, nobs, dOut, ndim);

产生的三维数据可用散布图显示：

% 画出 3 维资料的散布图 scatter3(x(1:nobs), x(nobs+1:2*nobs), x(2*nobs+1:3*nobs), ... 40.0, colours, 'filled'); % 加上文字、颜色等属性 dims = size(label); for i = 1 : nobs text(x(i), x(nobs+i), x(2*nobs+i), label(i,1:dims(2)), ... 'Color', [colours(i), colours(nobs+i), colours(2*nobs+i)], ... 'FontSize', 13); end

结果显示在图 3，图中可以观察出英国地区的数据 (黑点部分) 很接近蓝点 (terrestis 种)，再来才是红点 (sapidus 种)， 此意味着他们是属于同一类的。

图 3：田鼠资料度量尺度的散布图

本范例的 MATLAB 脚本位于 NAGToolboxDemos/Multivariate_methods/g03fa_demo.m， 也可下载 压缩文件。

G05 函数 - 随机数生成器

在许多的计算领域中，随机数序列的可靠性是一件很重要的任务，例如蒙地卡罗模拟。G05 函数包含了许多不同的算法，我们这里说明其中的两种算法。

在随机数生成的领域中有两个很明显的区别，一个为 伪随机 (pseudo)，一个为 准随机 (quasi)。 所谓的伪随机为统计属性上较趋近于 "真实" 的随机数 - 也就是说，这些随机数是由原始的随机数处理所取得 (计数器中所经过的时间频率等等)。 例如：在伪随机序列中每一个连续数字间的相关性是微乎其微的。相反的，准随机则没有这样的特性 - 他们会在空间中提供一个较为平均的分布值，这样的特性很适合应用于蒙地卡罗模拟中，我们给定序列的长度，会获得比伪随机更正确的估计值。

我们的例子可以更清楚的说明。以下的程序会使用 g05fa 函数产生两个伪随机序列：

% (x,y) 产生的数量 n = int32(10000); % 指定产生伪随机随机数 g05za('O'); % 在 [a, b] 区间中，产生长度为 n 的二个向量 a = 0; b = 1; [x] = g05fa(a, b, n); [y] = g05fa(a, b, n);

以下程序中采用 g05yd 函数的 Faure 算法，产生准随机的二维随机数序列：

% 指定随机数向量的维度 idim = int32(2); % 默认 Faure 算法 [iref, ifail] = g05yc(idim); % 产生二维向量 [quasi, irefOut, ifail] = g05yd(n, iref);

如果每个序列可以被解释为二维空间中的个别点，那么我们可用散布图表示 (如图 4)，其中可以很明显的显示出两种不同类型的随机数生成器的差异处。 虽然，他们似乎都可以随机的分布在二维的空间中，但是准随机序列以更一致的方式呈现 (技术上来说，这是因为产生序列中的每一点时，算法是以 最大化避免 的方式避免于其它数值接近)。

图 4：两种不同的随机数类别

本范例的 MATLAB 脚本位于 NAGToolboxDemos/Random_number_generation/g05fa_demo.m， 也可下载 压缩文件。