ACADEMIA

The NAG Toolbox makes the full functionality of the NAG Library available through MATLAB. An advantage of calling NAG via MATLAB is that many routine arguments become optional or unnecessary, which makes code easier to read and maintain. In addition, all of the NAG Library documentation (seventeen large volumes taking up a substantial amount of shelf space) has been converted to MATLAB help format, and so it is simple to access via MATLAB's usual documentation facilities. The documentation for each NAG Library routine includes example MATLAB code showing how to use the routine.

We have already highlighted some of the contents of the NAG Toolbox in our article here, we continue our exploration of the Toolbox, demonstrating how to call some of the other NAG routines, using MATLAB's plotting facilities to view the results.

Note: The code examples in this article have been extracted from demo scripts, and will not necessarily work properly if cut and pasted from this page into MATLAB. The full version of the scripts, which have been used to make the figures in this article, are available in this archive.

The C05 chapter - Finding the root of an equation

We start with the calculation of the zeros of a continuous function. The C05 chapter contains several routines for doing this (as well as other routines which solve a set of nonlinear equations in several unknowns). In our first example, we solve a quadratic equation.

  25x^2 - 10x + 1 = 0

Of course, it is trivial to determine analytically that both roots are at x = 0.2 (which provides a handy check on our numerical result - see below) but the C05 routines can be used to find the zeros of any function, as long as it is continuous. Here, we use routine c05ax; this evaluates the function via reverse communication, which allows us to display the intermediate results on the way to finding the root.

Here's the code, including some plotting commands:

% Specify the function as a MATLAB anonymous function - this facilitates
% changing the function for other examples.
myfun = @(x)25*x*x - 10*x + 1;

% Specify the initial estimate of the zero, and the value of the function
% there (the latter is only needed for plotting).
xstart = 1.0;
fstart = myfun(xstart);
xx = xstart;
fx = fstart;

% Specify other parameters for the NAG routine.  tol is the tolerance, 
% used by the routine to control the accuracy of the solution.
tol = 0.00001;
ir = int32(0);
c = zeros(26, 1);
ind = int32(1);

% Iterate until the solution is reached.
while (ind ~= int32(0))

% Call the NAG routine to get an improved estimate of the zero, 
% then plot it.
  [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 initial estimate, and final solution.
plot(axes, xstart, fstart, 'og', 'MarkerFaceColor', [0,1,0], 'MarkerSize', 8);
plot(axes, xx, fx, 'oy', 'MarkerFaceColor', [1,1,0], 'MarkerSize', 8);

and the results are shown in this picture:

As a second example, we can use the same code to find the root of a different function (for example, one that is harder to determine analytically) by changing the first few lines of the code above to:

% Specify the function as a MATLAB anonymous function - this facilitates
% changing the function for other examples.
myfun = @(x)x - exp(-x);

% Specify the initial estimate of the zero, and the value of the function
% there (the latter is only needed for plotting).
xstart = 5.0;

giving the results shown in Figure 2.

Figure 2: Finding the root of a transcendental function.

The MATLAB script for this demo is available as the file NAGToolboxDemos/Root_finding/c05ax_demo.m, distributed in this archive.

The E02 Chapter - Fitting a set of points with a cubic spline

A common task in many areas of science is to find a function which approximates a set of data points. The data may contain random perturbations - such as measurement errors - which it is desirable to smooth out. In our previous article, we looked at how to use the NAG Toolbox to compute a bicubic spline approximation to a two-dimensional set of data values given on a grid in the x-y plane; the result was displayed as a surface, f(x,y). Here, we investigate the one-dimensional case, where we seek to find a cubic spline curve f(x) which approximates a set of data values given on the x-axis. The NAG routine which does this is e02be, while the supplementary routine e02bb can be used to evaluate the spline at any point, given the results from e02be.

As in the two-dimensional example, a single parameter is specified to control the trade-off between closeness and smoothness of fit. A small value for this smoothing parameter, named s in e02be, will give a closer fit to the data, but the exact meaning of “small” in this context will vary according to the particular data set. In the extremes, a value of s = 0 will produce an interpolating spline, while a very large value will produce the weighted least-squares cubic polynomial approximation to the points.

% Here's the (x,y) data to be fitted, followed by the 
% weights for each point (set these to 1 if all points
% are equally important).  
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];

% Calculate work array lengths, according to the e02be documentation.
m = length(x);
nest = m+4;
lwrk = 4*m + 16*nest + 41;

% Now initialize some parameters required by the NAG routine.
start = 'C';
n = int32(0);
lamda = zeros(nest,1);
wrk = zeros(lwrk, 1);
iwrk = zeros(nest, 1, 'int32');

% Set the smoothness parameter and call the NAG routine to 
% compute the cubic spline approximation to the data.
s = 10000.0;
[nOut, lamdaOut, c, fp, wrkOut, iwrkOut, ifail] = ...
    e02be(start, x, y, w, s, n, lamda, wrk, iwrk);

% Set up the points at which the spline is to be evaluated.
px = [x(1) : 0.01 : x(m)];
pf = zeros(size(px));

% Now call the NAG routine e02bb to calculate the spline, and plot it.
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);

Figure 3 shows the data points and the spline for three values of s. It can be seen that, as s is decreased, the fit becomes less smooth, and closer to the points.

Figure 3: Fitting a cubic spline through data points.

The MATLAB script for this demo is available as the file NAGToolboxDemos/Curve_and_surface_fitting/e02be_demo.m, distributed in this archive.

The E01 Chapter - Interpolation through a set of points

In the previous example, when s = 0, we have a cubic spline curve that passes through all of the data points, and could be used to determine function values at any given x within the range of the supplied data points. This is interpolation, and the cubic spline can be a useful interpolating function for some data sets. However, as we have already seen in our discussion of the trade-off when determining values for the smoothness parameter, the interpolant can exhibit unwanted fluctuations in between the data points (see, for example, the Smoothness = 0 curve in Figure 3).

To illustrate this further, consider Figure 4, which shows a cubic spline interpolant, as fitted to another dataset.

Figure 4: Interpolating through data points using a piecewise cubic Hermite interpolant.

Here, we see that while the cubic spline appears to be a reasonable interpolant for low x values, it shows large fluctuations near the endpoints. These might be undesirable - for example, if the data represented some physical quantity whose value could not exceed 1, the behaviour of the cubic spline would be unacceptable.

However, other interpolants are available - for example, in chapter E01 of the NAG Toolbox, the e01be routine fits a piecewise cubic Hermite interpolant through a set of data points. This preserves monotonicity - that is, it will always increase, or decrease, over the same intervals as the data. As can be seen from Figure 4, this interpolant has no fluctuations, and might be more physically valid for our data set.

Here's the code for the e01be example. Once again, the NAG Toolbox provides an supplementary routine - here, e01bf - which can be used to evaluate the interpolant at any point, given the results from e01be.

% Here's the (x,y) data to be fitted.  Note that y is 
% monotonically increasing. 
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];
 
% Now call the NAG routine e01be to compute the Hermite interpolant
% through the data.
[d, ifail] = e01be(x, y);

% Set up the points at which the interpolant is to be calculated.
px = [x(1) : 0.01 : x(9)];

% Call the NAG routine e01bf to calculate the interpolant, and plot it.
[pf, ifail] = e01bf(x, f, d, px);
plot (px, pf, '-b', 'LineWidth', 2);

The MATLAB script for this demo is available as the file NAGToolboxDemos/Interpolation/e01be_demo.m, distributed in this archive.

The G05 Chapter - Creating a Gaussian copula

In this final section, we use the NAG Toolbox in an application that is somewhat more complicated. In statistics, a copula can be thought of as defining the correlation structure for a family of multivariate distributions. Each distribution in the family is constructed by "gluing" together two or more univariate distributions, with the copula supplying the "glue". Copulas can be applied to a wide range of simulation problems - for example, the so-called Gaussian copula is often used in financial modelling.

To illustrate, suppose we wish to simulate the joint distribution of several random variables. The kind of copula we use - Gaussian, Student's t, etc - describes the correlation structure between the variables, while a collection of univariate distributions (usually called the marginal distributions) define the distribution within each of the variables. Compared to more conventional multivariate distributions (such as, for example, the multivariate normal), the use of a copula allows each variable to have a different marginal distribution. Thus, for example, one variable may have a Gaussian (normal) marginal distribution, while another has a Beta distribution, a third is uniformly distributed, and so on.

In this example, we construct a bivariate distribution using a Gaussian copula with Beta marginal distributions. The input for the calculation of the multivariate distribution includes its covariance matrix, which contains details of the correlations between the different variables. The copula is created using NAG routine g05ra . The shape of the Beta distribution is defined by two parameters (usually called alpha and beta), and can range from the symmetric bell curve of the normal distribution (when alpha and beta are both 1) to being highly skewed or, at the other extreme, a uniform distribution. Here, we set alpha = 12.0 and beta = 5.0 for the first variable - which gives a skewed distribution - and alpha = beta = 5.0 for the second variable, which is distributed symmetrically. Here's the code:

% Prepare to call the NAG routine g05ra, which generates an array of pseudo-
% random variates from a Gaussian copula.  First, initialize some parameters 
% required by the routine.  
mode = int32(0);
igen = int32(1);
iseed = [int32(1); int32(2); int32(3); int32(4)];
r = zeros(7, 1);

% n is the number of values (actually value pairs) that will be generated.
n = int32(10000);

% Next, set the correlation coefficient for the two variables, and use it
% to construct the (upper triangle of the) covariance matrix of the
% distribution.  
corr = 0.8;
c = [1.0, corrl;
     0.0, 1.0;];

% Call g05ra to generate the array of variates.
[x, iseedOut, rOut, ifail] = g05ra(mode, c, n, igen, iseed, r);

% Now set the parameters for the two marginal distributions... 
alpha1 = 12.0;
beta1  = 5.0;
alpha2 = 5.0;
beta2  = 5.0;

% ...and call g01fe to compute the deviates from the Beta distribution.
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

The resultant distribution is then converted into a two-dimensional array of frequencies and displayed in MATLAB as a surface (though other representations - such as a contour plot - are also possible).

Figure 5: Creating a Gaussian copula.

It can be instructive to look at the shape of the copula for different marginal distributions, or for other correlations - for example, increasing the correlation coefficient causes the peak to become narrower. Furthermore, it's reasonably straightforward in MATLAB to construct an animation, based upon a loop over parameter values, which tracks the progress of the peak across the two-dimensional domain.

The MATLAB script for this demo is available as the file NAGToolboxDemos/Random_number_generation/g05ra_demo.m, distributed in this archive.