使用 NAG MATLAB® 工具箱

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

NAG MATLAB^®工具箱应用示例

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

NAG 工具箱能够让使用者在 MATLAB 的环境中，使用完整的 NAG 算法库函数。其中一个很大的特点是，透过 MATLAB 调用 NAG 函数，许多函数的参数都可以省略，或者以可选项方式呈现，这样可以让程序变的更容易阅读与维护。

所有 NAG 算法库的手册 - 几近 17 大册，如今已经全部转换为 MATLAB 的说明挡了。使用者可以很轻松的在 MATLAB 的环境中查找相关函数内容。在每一个函数中，也包含了 MATLAB 调用此函数的范例程序。

我们已经有其他文章说明如何的使用 NAG MATLAB 工具箱；我们持续在本文中探讨其他的函数，示范如何调用 NAG 函数，以及如何使用 MATLAB 的绘图函数显示计算结果。

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

C05 函数 - 方程式求根函数

我们先以求解连续函数的根为例。C05 包含了许多的函数求解，还有非线性方程的求解。在我们的第一个例子中，我们要求解一个二次方程式的根。

  25x^2 - 10x + 1 = 0

当然，我们可以透过分析判断出来，方程式的根是 x = 0.2 (以下，我们对结果进行验证)，C05 能够对任何的连续函数求根。这里，我们使用 c05ax 函数；它透过 reverse communication 方式求解，能够让我们显示求根的过程。

以下是部分的程序，包含了一些画图的命令：

% 设定 MATLAB 匿名函数
myfun = @(x)25*x*x - 10*x + 1;

% 设定根的默认估计值
xstart = 1.0;
fstart = myfun(xstart);
xx = xstart;
fx = fstart;

% 设定函数的其他参数。tol 变量为容许误差，提供函数控制求解的精确度
tol = 0.00001;
ir = int32(0);
c = zeros(26, 1);
ind = int32(1);

% 重复运行直到得到答案
while (ind ~= int32(0))

% 调用 NAG 函数取得逐渐逼近的根，并将值画出
  [xx, c, ind, ifail] = c05ax(xx, fx, tol, ir, c, ind);
  fx = myfun(xx);

  plot(axes, xx, fx, 'or', 'MarkerFaceColor', [1,0,0], 'MarkerSize', 8);

end

% 画出默认估计与最终结果
plot(axes, xstart, fstart, 'og', 'MarkerFaceColor', [0,1,0], 'MarkerSize', 8);
plot(axes, xx, fx, 'oy', 'MarkerFaceColor', [1,1,0], 'MarkerSize', 8);

将结果显示在图上：

图 1：求二次方程式的根

在第二个例子中，我们使用相同的程序求解不同的函数 (比较难直觉得计算出结果)，我们修改上面部分的程序片段：

% 设定 MATLAB 匿名函数
myfun = @(x)x - exp(-x);

% 设定根的默认估计值
xstart = 5.0;

结果显示在图 2

图 2：求二次方程式的根

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

在很多科学领域中，找出一个能够估计逼近所有数据点的函数是极为普遍的问题。这些数据也许包含许多混乱随机的数据，例如量测的误差等，我们便需要将其平滑处理的。在先前的文章中，我们看到了如何使用 NAG 工具箱对给定的二维 x-y 网格进行 bicubic spline 的估计，并将其结果以 f(x,y) 的曲面图表示。这里，我们来探讨一维的例子，我们来求取平面上 f(x) 的曲线拟合。NAG 的 e02be 函数可以求解这样的问题，与其搭配的函数 e02bb 可依据 e02be 函数的结果求取任何数值的拟合结果。

在二维的例子中，我们用单一的参数来设定拟合紧密与平滑的程度。在 e02be 函数中，采用 s 参数来设定拟合的程度，然而但是其确切意义会随着数据特性的不同而有所改变。在极端的例子下，当s = 0 的参数设定时会形成 插值样条；若参数值设的很大，会对数据点产生 加权最小二乘的三次多项式 估计。

% 指定 (x,y) 要拟合的数据，并对每一数据点指定其比重
% 如果所有的数据点比重都一样的话，可以都设为 1
x = [0.5; 1; 1.5; 2; 2.5; 3; 4; 4.5; 5; 5.5; 6; 7; 7.5];
y = [-0.372; 0.431; 1.69; 2.11; 3.1; 4.23; 4.35;
     4.81; 4.61; 4.79; 5.23; 6.35; 7.19];
w = [2; 1.5; 1; 3; 1; 0.5; 1; 2; 2.5; 1; 3; 1; 2];

% 计算工作矩阵的大小，参考 e02be 文件说明
m = length(x);
nest = m+4;
lwrk = 4*m + 16*nest + 41;

% 设定 NAG 函数的相关参数
start = 'C';
n = int32(0);
lamda = zeros(nest,1);
wrk = zeros(lwrk, 1);
iwrk = zeros(nest, 1, 'int32');

% 设定平滑参数并调用 NAG 函数，计算 cubic spline 估计
s = 10000.0;
[nOut, lamdaOut, c, fp, wrkOut, iwrkOut, ifail] = ...
    e02be(start, x, y, w, s, n, lamda, wrk, iwrk);

% 设定要进行 spline 求值的数据
px = [x(1) : 0.01 : x(m)];
pf = zeros(size(px));

% 调用 NAG e02bb 函数，计算其 spline 值，并将结果画出
for i = 1 : length(px)
  [pf(i), ifail] = e02bb(lamdaOut, c, px(i), 'ncap7', nOut);
end
plot(px, pf, 'Color', [.8 0 0], 'Linewidth', 2);

图 3 显示数据点及不同 s 参数设定后的拟合图形。您可以发现，当 s 值减少时，拟合的曲线就会越不平滑，而且会接近数据点。

图 3：利用 cubic spline 拟合数据点

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

E01 函数 - 数据点的插值函数

在前面的例子中，当 s = 0 时，我们会得到通过每一个数据点的 cubic spline 曲线，而且能够用来决定在数据给定的区间中任何 x 值的函数值。这便叫做插值，对许多资料而言，cubic spline 是一个很有用的插值函数。无论如何，在我们先前的讨论中我们可以了解到平滑参数设定的相关得失，插值可以表现数据点间的波动。(例如：在图 3 中设定 Smoothness = 0 的曲线)。

为了更进一步的说明，在图 4 中我们对另一组数据进行 cubic spline 插值。

图 4：使用分段三次 Hermite 插值法。

此例中，在 x 的值较小的地方，我们可以看到 cubic spline 插值表现的相当理想，在末端的地方，确显现出很多巨大的波动。这样的结果可能很不理想，举个例来说，若某些实际数据不可能超过 1，那么 cubic spline 所呈现的结果就无法被我们所接受了。

然而，我们也提供了其他不同的插值函数 - 例如 E01 中的 e01be 函数，可以对数据点进行分段三次插补 (Piecewise Cubic Hermite Interpolant)。它确保插值的单一性的特性，也就是说它计算的结果总是会在数据的区间中递增或递减。我们可以从图 4 中看到这样的结果，这样的插值结果不会形成波动，也更符合我们数据的特性。

以下是 e01be 的程序片段。同样的，NAG 工具箱也提供补充函数，e01bf 函数能够依据 e01be 函数的结果，计算每一个指定点的內插。

% (x,y) 数据，y 的值是递增的
x = [7.5;
     8.09;
     8.19;
     8.69;
     9.19;
     10;
     12;
     15;
     20];
y = [0;
     0;
     0.04375;
     0.16918;
     0.46943;
     0.94374;
     0.99864;
     0.99992;
     0.99999];

% 调用 NAG e01be 函数，计算 Hermite 插值
[d, ifail] = e01be(x, y);

% 设定要求插值的数据点
px = [x(1) : 0.01 : x(9)];

% 调用 NAG e01bf 函数，计算插值并绘图
[pf, ifail] = e01bf(x, f, d, px);
plot (px, pf, '-b', 'LineWidth', 2);

本范例的 MATLAB 脚本位于 NAGToolboxDemos/Interpolation/e01be_demo.m，可以下载压缩文件。

G05 函数 - 计算 Gaussian copula

最后一部分，我们将介绍一个较为复杂的 NAG 工具箱函数的应用。在统计应用中，copula 会被应用在多变量分布的相关性结构定义中。在这个问题中的每一个分布都是由二个或数个的单变数的分布所结合的，而 copula 便提供这样的功能。Copulas 很普遍的应用在模拟问题上。例如：一般称的 Gaussian copula 被广泛的应用在金融模型建立中。

为了详细说明，我们假设我们要模拟几个随机数的共同分布。我们使用 Gaussian、Student's t 的 Copula，来描述这些变量间的相关性结构，集合每一个变量间的单变量分布 (通常称为 边界分布)。与更正规花的多变量分布 (例如：多变量常态分配) 比较，透过 copula 的使用，可以让每一个变量有不同的边际分布。例如：一个变量是 Gaussian (常态) 边际分布，同时另一个变量是 Beta 分布；第三个变数是 uniformly 分布等等。

在我们的例子中，我们使用 Gaussian coupla 与 Beta 边际分布建构双变量的分布。多变量分布的计算中需要输入共变异矩阵，它包含了各个不同变量的相关系数。我们使用 NAG g05ra 函数计算 copula。 Beta 分布的形状是由两个参数所定义的 (通常是称 alpha 与 beta)，当 alpha 与 beta 都是 1 时，他们是对称的钟形曲线，或者是相当偏斜的形状，或者是单峯分布。我们在这里对第一个变量设定 alpha = 12.0，beta = 5.0 - 这会形成偏态分布；另外对第二个变量设定 alpha = beta = 5.0，是一个对称的分布。以下是我们的程序：

% 调用 NAG g05ra 函数预备工作，从 Gaussian copula 产生随机数
% 默认化函数调用所需的参数
mode = int32(0);
igen = int32(1);
iseed = [int32(1); int32(2); int32(3); int32(4)];
r = zeros(7, 1);

% n 是所需要产生的数量
n = int32(10000);

% 为两个变量设定其相关系数，并依此建构其方差矩阵 (上三角)
corr = 0.8;
c = [1.0, corrl;
     0.0, 1.0;];

% 调用 g05ra 产生变量
[x, iseedOut, rOut, ifail] = g05ra(mode, c, n, igen, iseed, r);

% 设定两个边际分布的参数...
alpha1 = 12.0;
beta1  = 5.0;
alpha2 = 5.0;
beta2  = 5.0;

% ...调用 g01fe 计算其 Beta 分布的偏差值
tol = 1.0;
for ii = 1 : n
  [x(ii,1)] = g01fe(x(ii,1), alpha1, beta1, tol);
  [x(ii,2)] = g01fe(x(ii,2), alpha2, beta2, tol);
end

由此产生的分布转换为二维的频率矩阵，并以 MATLAB 的 surface 绘图指令画出。

图 5：Gaussian copula

这个例子可以启发您让您了解不同边际分布或相关系数的 copula 型态 - 例如：增加相关系数会导致顶端更为狭窄。此外，可以在 MATLAB 中，透过循环中参数的调整，很直觉的建构出一个动画图形，观察顶端在二维空间中的变化状况。

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