C# Weighted Least Squares Example

using System;

using CenterSpace.NMath.Core;


namespace CenterSpace.NMath.Examples.CSharp
{
  /// <summary>
  /// A .NET example in C# demonstrating the features of the classes for solving weighted least
  /// squares problems.
  /// </summary>
  class WeightedLeastSquaresExample
  {

    static void Main( string[] args )
    {
      // Construct some test data.
      var A = new DoubleMatrix( "5x2[1 2  1 3  1 6  1 10  1 7]" );
      var b = new DoubleVector( "[ 3 6 8 10 11]" );

      // Constructs some arbitrary weights. Observations are given 
      // increasingly higher weights.
      var weights = new DoubleVector( A.Rows, .2, .2 );

      Console.WriteLine();

      Console.WriteLine( "Solving weighted least squares problem, Ax = b " );
      Console.WriteLine();
      Console.WriteLine( "with A = " );
      Console.WriteLine( A.ToTabDelimited() );
      Console.WriteLine( "and b = " + b );
      Console.WriteLine();
      Console.WriteLine( "with weights = " + weights );
      Console.WriteLine();

      // Construct a weighted least squares object and solve the problem. The third
      // parameter says do not prepend a column of ones to the input matrix A (the
      // column of ones represents a constant term in the model) as the first column
      // of A is already ones.
      var wls = new DoubleCOWeightedLeastSq( A, weights, false );
      DoubleVector solution = wls.Solve( b );
      Console.WriteLine( "Solution = " + solution.ToString( "G5" ) );
      Console.WriteLine( "Residuals = " + wls.ResidualVector( b ).ToString( "G5" ) );

      // Construct a random problem, obtain an ordinary least squares solution,
      // then use the resulting residuals to solve the same problem using a
      // weighted least squares - downweighting the outliers.
      var rng = new RandGenUniform( -1, 2, 0x24 );
      int rows = 10;
      int cols = 3;
      A = new DoubleMatrix( rows, cols, rng );
      b = new DoubleVector( rows, rng );

      // Ordinary Least Squares (OLS):
      var ols = new DoubleQRLeastSq();
      ols.Factor( A );
      DoubleVector olsSolution = ols.Solve( b );
      DoubleVector residuals = ols.ResidualVector( b );
      Console.WriteLine();
      Console.WriteLine( "OLS solution = " + olsSolution.ToString( "G5" ) );

      // Now, construct a Bisquare weighting function object
      // and use it to compute weights:
      var weightingFunction = new DoubleBisquareWeightingFunction( A );
      weightingFunction.GetWeights( residuals, ref weights );

      // Perform the factorization with the input matrix and weights. Once the factorization is
      // done we can use the weighted least squares to quickly solve for different right hand
      // sides. The third parameter to the Factor method indicates that we want a column of 
      // ones prepended to the input matrix A representing a constant term in the model.
      wls.Factor( A, weights, true );

      // Now solve the weighted least squares problem.
      DoubleVector wlsSolution1 = wls.Solve( b );
      Console.WriteLine( "WLS solution = " + wlsSolution1.ToString( "G5" ) );

      // Solve with a different right hand side.
      var c = new DoubleVector( rows, rng );
      DoubleVector wlsSolution2 = wls.Solve( c );
      Console.WriteLine( "WLS solution 2 = " + wlsSolution2.ToString( "G5" ) );

      Console.WriteLine();
      Console.WriteLine( "Press Enter Key" );
      Console.Read();
    }
  }
}