Example: Multiquadric Radial Basis Function Approximation

Data is generated from the function
e^\frac{y}{2.0}\sin{x}\cos\frac{y}{2.0}
where a number of (x,y ) pairs make up a set of randomly chosen points. Random noise is added to the values, a Hardy multiquadric radial basis function is specified \sqrt{r^2+\delta^2} and a radial basis approximation of the noisy data is computed. The radial basis fit is then compared to the original function at another set of randomly chosen points. Both the average normalized error and the maximum normalized error are computed and printed.

In this example, the parameter of the Hardy multiquadric radial basis function \delta = 5.5. The function is sampled at 100 random points and the error is computed at 10000 random points.

using System;
using Imsl.Math;

public class RadialBasisEx3
{
    
    public static void  Main(String[] args)
    {
        int nDim = 2;
        
        // Sample, with noise, the function at 100 randomly chosen points
        int nData = 100;
        double[,] xData = new double[nData, nDim];
        double[] fData = new double[nData];
        double[] row = new double[nDim];

        Imsl.Stat.Random rand = new Imsl.Stat.Random(123457);
        rand.Multiplier = 16807;
        double[] noise = new double[nData * nDim];
        for (int k = 0; k < nData; k++)
        {
            for (int i = 0; i < nDim; i++)
            {
                noise[k * 2 + i] = 1.0d - 2.0d * (double) rand.NextDouble();
                xData[k, i] = 3 * noise[k * 2 + i];
            }
            // noisy sample
            for(int j = 0; j<nDim; j++)
                row[j]=xData[k,j];
            fData[k] = fcn(row) + noise[k * 2] / 10;
        }
        
        // Compute the radial basis approximation using 100 centers
        int nCenters = 100;
        RadialBasis rb = new RadialBasis(nDim, nCenters);
        rb.RadialFunction = new RadialBasis.HardyMultiquadric(5.5);
        rb.Update(xData, fData);
        
        // Compute the error at a randomly selected set of points
        int nTest = 10000;
        double maxError = 0.0;
        double aveError = 0.0;
        double maxMagnitude = 0.0;
        Imsl.WarningObject w = Imsl.Warning.WarningObject;
        Imsl.Warning.WarningObject = null;
        double[][] x = new double[nTest][];
        for (int i2 = 0; i2 < nTest; i2++)
        {
            x[i2] = new double[nDim];
        }
        noise = new double[nTest * nDim];
        
        for (int i = 0; i < nTest; i++)
        {
            for (int j = 0; j < nDim; j++)
            {
                noise[i * 2 + j] = 1.0d - 2.0d * rand.NextDouble();
                x[i][j] = 3 * noise[i * 2 + j];
            }
            double error = Math.Abs(fcn(x[i]) - rb.Eval(x[i]));
            maxMagnitude = Math.Max(Math.Abs(fcn(x[i])),
                maxMagnitude);
            aveError += error;
            maxError = Math.Max(error, maxError);
        }
        aveError /= nTest;
        
        Imsl.Warning.WarningObject = w;
        Console.WriteLine("Average normalized error is " +
            aveError / maxMagnitude);
        Console.WriteLine("Maximum normalized error is " +
            maxError / maxMagnitude);
    }
    
    // The function to approximate
    internal static double fcn(double[] x)
    {
        return Math.Exp((x[1]) / 2.0) * Math.Sin(x[0]) -
            Math.Cos((x[1]) / 2.0);
    }
}

Output

Average normalized error is 0.0110861431448538
Maximum normalized error is 0.054260728660165
Link to C# source.