Example: "Custom" 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 "custom" Polyharmonic Spline radial basis function is specified
\varphi(r)=\left\{\begin{array}{ll}r^{k}&k=1,3,5, \ldots\\r^{k}\ln{r}&k=2,4,6,\ldots\end{array}\right.
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 order of the Polyharmonic Spline, k =2. The function is sampled at 200 random points and the error is computed at 10000 random points.

using System;
using Imsl.Math;

public class RadialBasisEx2
{
    
    // The function to approximate
    static double fcn(double[] x)
    {
        return Math.Exp((x[1]) / 2.0) * Math.Sin(x[0]) -
            Math.Cos((x[1]) / 2.0);
    }
    
    public class PolyHarmonicSpline : RadialBasis.IFunction
    {
        virtual public int Order
        {
            get
            {
                return order;
            }
            
        }
        virtual public bool EvenOrder
        {
            get
            {
                return isEven;
            }
            
        }
        
        private int order = 3;
        private bool isEven = false;
        
        public PolyHarmonicSpline(int order)
        {
            this.isEven = order % 2 == 0;
            this.order = order;
        }
        
        public virtual double F(double x)
        {
            if (this.isEven)
            {
                return Math.Pow(x, order) * Math.Log(x);
            }
            return Math.Pow(x, order);
        }
        
        public virtual double G(double x)
        {
            if (order == 1)
            {
                return 1;
            }
            if (this.isEven)
            {
                return order * Math.Pow(x, order - 1) * Math.Log(x) +
                    Math.Pow(x, order - 1);
            }
            return order * Math.Pow(x, order - 1);
        }
    }


    public static void  Main(String[] args)
    {
        int nDim = 2;
        
        // Sample, with noise, the function at 100 randomly chosen points
        int nData = 200;
        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 * 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 PolyHarmonicSpline(2);
        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;
        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 * (double) 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;
        
        Console.WriteLine("Average normalized error is " + aveError /
            maxMagnitude);
        Console.WriteLine("Maximum normalized error is " + maxError /
            maxMagnitude);
        Console.WriteLine("Using even order equation: " +
            ((PolyHarmonicSpline) rb.RadialFunction).EvenOrder);
    }
}

Output

Average normalized error is 0.0180558727060151
Maximum normalized error is 0.257565310407464
Using even order equation: True
Link to C# source.