C# Weighted Least Squares Example

using System;

using CenterSpace.NMath.Core;
using CenterSpace.NMath.Matrix;

namespace CenterSpace.NMath.Matrix.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.
      DoubleMatrix A = new DoubleMatrix("5x2[1 2  1 3  1 6  1 10  1 7]");
      DoubleVector b = new DoubleVector("[ 3 6 8 10 11]");

      // Constructs some arbitrary weights. Observations are given 
      // increasingly higher weights.
      DoubleVector weights = new DoubleVector(A.Rows, .2, .2);
      Console.WriteLine("Solving weighted least squares problem, Ax = b ");
      Console.WriteLine("with A = ");
      Console.WriteLine(A.ToTabDelimited());
      Console.WriteLine("and b = " + b);
      Console.WriteLine("With weights = " + weights);

      // 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.
      DoubleCOWeightedLeastSq wls = new DoubleCOWeightedLeastSq(A, weights, false);
      DoubleVector solution = wls.Solve(b);
      Console.WriteLine("Solution = " + solution);
      Console.WriteLine("Residuals = " + wls.ResidualVector(b));

      // 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.
      RandGenUniform 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):
      DoubleQRLeastSq ols = new DoubleQRLeastSq();
      ols.Factor(A);
      DoubleVector olsSolution = ols.Solve(b);
      DoubleVector residuals = ols.ResidualVector(b);
      Console.WriteLine();
      Console.WriteLine("OLS solution = " + olsSolution);

      // Now, construct a Bisquare weighting function object
      // and use it to compute weights:
      DoubleBisquareWeightingFunction 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);

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

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