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

NAG MATLAB®工具箱应用示例

使用 NAG MATLAB® 工具箱 - Part 3

本系列的文章可参考 http://tw.nag-gc.com/IndustryArticles/usingtoolboxmatlab.asp 以及 http://tw.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， 也可下载 压缩文件。