Blog

Archive for the ‘NMath Tutorial’ Category

NMath Tutorial videos

Thursday, October 31st, 2013

We are proud to announce a series of NMath tutorial videos on how to use CenterSpace’s math library products. We are starting, naturally, with Getting Started with NMath.

You can download it here: MP4

Please let us know which topics you want us to cover. Email support@centerspace.net

Cheers,
Trevor

Distribution Fitting Demo

Monday, April 9th, 2012

A customer recently asked how to fit a normal (Gaussian) distribution to a vector of experimental data. Here’s a demonstration of how to do it.

Let’s start by creating a data set: 100 values drawn from a normal distribution with known parameters (mean = 0.5, variance = 2.0).

int n = 100;
double mean = .5;
double variance = 2.0;
var data = new DoubleVector( n, new RandGenNormal( mean, variance ) );

Now, compute y values based on the empirical cumulative distribution function (CDF), which returns the probability that a random variable X will have a value less than or equal to x–that is, f(x) = P(X <= x). Here’s an easy way to do, although not necessarily the most efficient for larger data sets:

var cdfY = new DoubleVector( data.Length );
var sorted = NMathFunctions.Sort( data );
for ( int i = 0; i < data.Length; i++ )
{
  int j = 0;
  while ( j < sorted.Length && sorted[j] <= data[i] ) j++;
  cdfY[i] = j / (double)data.Length;
}

The data is sorted, then for each value x in the data, we iterate through the sorted vector looking for the first value that is greater than x.

We’ll use one of NMath’s non-linear least squares minimization routines to fit a normal distribution CDF() function to our empirical CDF. NMath provides classes for fitting generalized one variable functions to a set of points. In the space of the function parameters, beginning at a specified starting point, these classes finds a minimum (possibly local) in the sum of the squared residuals with respect to a set of data points.

A one variable function takes a single double x, and returns a double y:

y = f(x)

A generalized one variable function additionally takes a set of parameters, p, which may appear in the function expression in arbitrary ways:

y = f(p1, p2,..., pn; x)

For example, this code computes y=a*sin(b*x + c):

public double MyGeneralizedFunction( DoubleVector p, double x )
{
  return p[0] * Math.Sin( p[1] * x + p[2] );
}

In the distribution fitting example, we want to define a parameterized function delegate that returns CDF(x) for the distribution described by the given parameters (mean, variance):

Func<DoubleVector, double, double> f =
  ( DoubleVector p, double x ) =>
    new NormalDistribution( p[0], p[1] ).CDF( x );

Now that we have our data and the function we want to fit, we can apply the curve fitting routine. We’ll use a bounded function fitter, because the variance of the fitted normal distribution must be constrained to be greater than 0.

var fitter = new BoundedOneVariableFunctionFitter<TrustRegionMinimizer>( f );
var start = new DoubleVector( new double[] { 0.1, 0.1 } );
var lowerBounds = new DoubleVector( new double[] { Double.MinValue, 0 } );
var upperBounds = 
   new DoubleVector( new double[] { Double.MaxValue, Double.MaxValue } );
var solution = fitter.Fit( data, cdfY, start, lowerBounds, upperBounds );
var fit = new NormalDistribution( solution[0], solution[1] );
 
Console.WriteLine( "Fitted distribution:\nmean={0}\nvariance={1}",
  fit.Mean, fit.Variance );

The output for one run is

Fitted distribution: 
mean=0.567334190790594
variance=2.0361207956132

which is a reasonable approximation to the original distribution (given 100 points).

We can also visually inspect the fit by plotting the original data and the CDF() function of the fitted distribution.

ToChart( data, cdfY, SeriesChartType.Point, fit,
  NMathStatsChart.DistributionFunction.CDF );
 
private static void ToChart( DoubleVector x, DoubleVector y,
  SeriesChartType dataChartType, NormalDistribution dist,
  NMathStatsChart.DistributionFunction distFunction )
{
  var chart = NMathStatsChart.ToChart( dist, distFunction );
  chart.Series[0].Name = "Fit";
 
  var series = new Series() {
    Name = "Data",
    ChartType = dataChartType
  };
  series.Points.DataBindXY( x, y );
  chart.Series.Insert( 0, series );
 
  chart.Legends.Add( new Legend() );
  NMathChart.Show( chart );
}

CDF() of fitted distribution

We can also look at the probability density function (PDF) of the fitted distribution, but to do so we must first construct an empirical PDF using a histogram. The x-values are the midpoints of the histogram bins, and the y-values are the histogram counts converted to probabilities, scaled to integrate to 1.

int numBins = 10;
var hist = new Histogram( numBins, data );
 
var pdfX = new DoubleVector( hist.NumBins );
var pdfY = new DoubleVector( hist.NumBins );
for ( int i = 0; i < hist.NumBins; i++ )
{
  // use bin midpoint for x value
  Interval bin = hist.Bins[i];
  pdfX[i] = ( bin.Min + bin.Max ) / 2;
 
   // convert histogram count to probability for y value
   double binWidth = bin.Max - bin.Min;
   pdfY[i] = hist.Count( i ) / ( data.Length * binWidth );
}
 
ToChart( pdfX, pdfY, SeriesChartType.Column, fit,
  NMathStatsChart.DistributionFunction.PDF );

PDF() of fitted distribution

You might be tempted to try to fit a distribution PDF() function directly to the histogram data, rather than using the CDF() function like we did above, but this is problematic for several reasons. The bin counts have different variability than the original data. They also have a fixed sum, so they are not independent measurements. Also, for continuous data, fitting a model based on aggregated histogram counts, rather than the original data, throws away information.

Ken

Clearing a vector

Wednesday, November 9th, 2011

A customer recently asked us for the best method to zero out a vector. We decided to run some tests to find out. Here are the five methods we tried followed by performance timing and any drawbacks.

The following tests were performed on a DoubleVector of length 10,000,000.

1) Create a new vector. This isn’t really clearing out an existing vector but we thought we should include it for completeness.

1
 DoubleVector v2 = new DoubleVector( v.Length, 0.0 );

The big drawback here is that you’re creating new memory. Time: 419.5ms

2) Probably the first thing to come to mind is to simply iterate through the vector and set everything to zero.

1
2
3
4
for ( int i = 0; i < v.Length; i++ )
{
  v[i] = 0.0;
}

We have to do some checking in the index operator. No new memory is created. Time: 578.5ms

3) In some cases, you could iterate through the underlying array of data inside the DoubleVector.

1
2
3
4
for ( int i = 0; i &lt; v.DataBlock.Data.Length; i++ )
{
  v.DataBlock.Data[i] = 0.0;
}

This is a little less intuitive. And, very importantly, it will not work with many views into other data structures. For example, a row slice of a matrix. However, it’s easier for the CLR to optimize this loop. Time: 173.5ms

4) We can use the power of Intel’s MKL to multiply the vector by zero.

1
 v.Scale( 0.0 );

Scale() does this in-place. No new memory is created. In this example, we assume that MKL has already been loaded and is ready to go which is true if another MKL-based NMath call was already made or if NMath was initialized. This method will work on all views of other data structures. Time: 170ms

5) This surprised us a bit but the best method we could find was to clear out the underlying array using Array.Clear() in .NET

1
 Array.Clear( v.DataBlock.Data, 0, v.DataBlock.Data.Length );

This creates no new memory. However, this will not work with non-contiguous views. However, this method is very fast. Time: 85.8ms

To make efficiently clearing a vector simpler for NMath users we have created a Clear() method and a Clear( Slice ) method on the vector and matrix classes. It will do the right thing in the right circumstance. It will be released in NMath 5.2 in 2012.

Test Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
using System;
using CenterSpace.NMath.Core;
 
namespace Test
{
  class ClearVector
  {
    static int size = 100000000;
    static int runs = 10;
    static int methods = 5;
 
    static void Main( string[] args )
    {
      System.Diagnostics.Stopwatch sw = new System.Diagnostics.Stopwatch();
      DoubleMatrix times = new DoubleMatrix( runs, methods );
      NMathKernel.Init();
 
      for ( int run = 0; run < runs; run++ )
      {
        Console.WriteLine( "Run {0}...", run );
        DoubleVector v = null;
 
        // Create a new one
        v = new DoubleVector( size, 1.0, 2.0 );
        sw.Start();
        DoubleVector v2 = new DoubleVector( v.Length, 0.0 );
        sw.Stop();
        times[run, 0] = sw.ElapsedMilliseconds;
        Console.WriteLine( Assert( v2 ) );
 
        // iterate through vector
        v = new DoubleVector( size, 1.0, 2.0 );
        sw.Reset();
        sw.Start();
        for ( int i = 0; i < v.Length; i++ )
        {
          v[i] = 0.0;
        }
        sw.Stop();
        times[run, 1] = sw.ElapsedMilliseconds;
        Console.WriteLine( Assert( v ) );
 
        // iterate through array
        v = new DoubleVector( size, 1.0, 2.0 );
        sw.Reset();
        sw.Start();
        for ( int i = 0; i < v.DataBlock.Data.Length; i++ )
        {
          v.DataBlock.Data[i] = 0.0;
        }
        sw.Stop();
        times[run, 2] = sw.ElapsedMilliseconds;
        Console.WriteLine( Assert( v ) );
 
        // scale
        v = new DoubleVector( size, 1.0, 2.0 );
        sw.Reset();
        sw.Start();
        v.Scale( 0.0 );
        sw.Stop();
        times[run, 3] = sw.ElapsedMilliseconds;
        Console.WriteLine( Assert( v ) );
 
        // Array Clear
        v = new DoubleVector( size, 1.0, 2.0 );
        sw.Reset();
        sw.Start();
        Array.Clear( v.DataBlock.Data, 0, v.DataBlock.Data.Length );
        sw.Stop();
        times[run, 4] = sw.ElapsedMilliseconds;
        Console.WriteLine( Assert( v ) );
        Console.WriteLine( times.Row( run ) );
      }
      Console.WriteLine( "Means: " + NMathFunctions.Mean( times ) );
    }
 
    private static bool Assert( DoubleVector v )
    {
      if ( v.Length != size )
      {
        return false;
      }
      for ( int i = 0; i < v.Length; ++i )
      {
        if ( v[i] != 0.0 )
        {
          return false;
        }
      }
      return true;
    }
  }
}

– Trevor

Initializing NMath

Wednesday, November 9th, 2011

NMath uses Intel’s Math Kernel Library (MKL) internally. This code contains native, optimized code to wring out the best performance possible.

There is a one-time delay when the appropriate x86 or x64 native code is loaded. This cost can be easily controlled by the developer by using the NMathKernel.Init() method. Please see Initializing NMath for more details.

– Trevor

Fitting Geometric Primitives to Points Using Nonlinear Least Squares

Tuesday, August 9th, 2011

We were recently contacted by a customer looking for help on how to use NMath to fit geometric primitives to clouds of 2D points. The solution is to cast the problem as a minimization. In the space of the parameters which define the geometric object, minimize the residuals with respect to the 2D points. NMath provides .NET classes for solving nonlinear least squares problems such as this, using the Trust-Region or Levenberg-Marquardt methods.

For instance, let’s try fitting a circle to a set of points. A circle is defined by three parameters: a center (a,b), and a radius r. The circle is all points such that (x-a)^2 + (y – b)^2 = r^2. To setup the minimization problem, we first need to define a function which given a set of circle parameters, returns the residuals with respect to a cloud of 2D points. There are several methods for doing this in NMath. In the following C# code, we extend the abstract base class DoubleMultiVariableFunction, and override the Evaluate() method:

public class CircleFitFunction : DoubleMultiVariableFunction
{
 
  public CircleFitFunction( DoubleVector x, DoubleVector y )
    : base( 3, x.Length )
  {
    if( x.Length != y.Length )
      throw new Exception( "Unequal number of x,y values." );
 
    X = x;
    Y = y;
  }
 
  public DoubleVector X { get; internal set; }
  public DoubleVector Y { get; internal set; }
 
  public override void Evaluate( DoubleVector parameters, ref DoubleVector residuals )
  {
    // parameters of circle with center (a,b) and radius r
    double a = parameters[0];
    double b = parameters[1];
    double r = parameters[2];
 
    for( int i = 0; i < X.Length; i++ )
    {
      // distance of point from circle center
      double d = Math.Sqrt( Math.Pow( X[i] - a, 2.0 ) + Math.Pow( Y[i] - b, 2.0 ) );
 
      residuals[i] = d - r;
    }
  }
}

The constructor takes the set of x,y points to fit, which are stored on the class. Note that we call the base constructor with the dimensions of the domain and range of the function:

: base( 3, x.Length )

In this case, the function maps the 3-dimensional space of the circle parameters to the x.Length-dimension space of the residuals.

Next we override the Evaluate() method, which takes an input vector and output vector passed by reference. Our implementation computes the residual for each point with respect to a circle defined by the given parameters. We calculate the distance, d, of each point from the circle center; the residual is then equal to d – r. The nonlinear least squares method will minimize the L2 norm of this function.

To demonstrate the fitting process, let’s first start with a set of points generated from a circle of known parameters. For example, this C# code creates 20 x,y points evenly spaced around a circle with center (0.5, 0.25) and radius 2, plus some added noise:

int n = 20;
DoubleVector x = new DoubleVector( n );
DoubleVector y = new DoubleVector( n );
 
double a = 0.5;
double b = 0.25;
double r = 2;
 
RandGenUniform noise = new RandGenUniform( -.75, .75 );
for( int i = 0; i < n; i++ )
{
  double t = i * 2 * Math.PI / n;
  x[i] = a + r * Math.Cos( t ) + noise.Next();
  y[i] = b + r * Math.Sin( t ) + noise.Next();
}

Now that we have our function defined and some points to fit, performing the minimization takes only a few lines of code. This C# code finds the minimum and prints the solution and final residual:

CircleFitFunction f = new CircleFitFunction( x, y );
TrustRegionMinimizer minimizer = new TrustRegionMinimizer();
DoubleVector start = new DoubleVector( "0.1 0.1 0.1" );
DoubleVector solution = minimizer.Minimize( f, start );
 
Console.WriteLine( "solution = " + solution );
Console.WriteLine( "final residual = " + minimizer.FinalResidual );

Sample output:

solution = [ 0.446122611523468 0.175483486509012 1.93511389538286 ]
final residual = 1.62414259527024

Starting from point (0.1, 0.1, 0.1) in the parameter space of the circle, we minimized the sum of the squared residuals with respect to the (x,y) points. In the run above, the minimum found was center (0.45, 0.18) and radius 1.9, close to the actual parameters which generated our (noisy) data. To visually inspect the fit, we can use NMath with the Microsoft Chart Controls for .NET.

Fitted Circle

Now that we’re sure the minimization is working as desired, let’s try it again with some random points:

int n = 10;
RandGenUniform rng = new RandGenUniform( -1, 2 );
DoubleVector x = new DoubleVector( n, rng );
DoubleVector y = new DoubleVector( n, rng );
 
CircleFitFunction f = new CircleFitFunction( x, y );
TrustRegionMinimizer minimizer = new TrustRegionMinimizer();
DoubleVector start = new DoubleVector( "0.1 0.1 0.1" );
DoubleVector solution = minimizer.Minimize( f, start );

Again, plotting the solution:

Fitted CircleIn some runs, the fit may seem counter-intuitive. For example:

Fitted Circle

But if you think about it, this sort of fit makes perfect sense, given how we defined our fit function. We’re asking the minimizer to minimize the residuals with respect to the circle perimeter, and allowing it vary the circle center and radius however it can to achieve that goal. In an extreme case, if our points fall approximately along a line, the best fit will be a circle with a very large radius, so that the circle perimeter is nearly linear as it passes through the points.

Fitted CircleIn this case, the fitted circle has center (428.0, -757.8) and radius 870.6.  If we wish to avoid solutions such as these, we could define our fit function differently. For example, we might constrain the circle center to be the center of mass of the points, and vary only the radius.

Obviously, a similar technique can be used to fit other geometric primitives, though as the complexity of the shape increases, so does the complexity of the fit function.

Ken