|
///////////////////////////////////////////////////////////////////////////////
//
// Program.cs
//
// By Philip R. Braica (HoshiKata@aol.com, VeryMadSci@gmail.com)
//
// Distributed under the The Code Project Open License (CPOL)
// http://www.codeproject.com/info/cpol10.aspx
///////////////////////////////////////////////////////////////////////////////
// Using.
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
// Namespace.
namespace SolverDemo
{
/// <summary>
/// Linear least squares for 2 variables, 1st and second order,
/// beyond that requires a decent matrix math package with good pseudo-inverse.
/// Note that non-linear is more accurate when there is noise but overkill for
/// this application, this is just a tiny quick fitter for use by the Solver class.
///
/// See the following for more information:
/// http://en.wikipedia.org/wiki/Simple_linear_regression
/// http://en.wikipedia.org/wiki/Invertible_matrix
/// </summary>
public class LinearLeastSquares
{
/// <summary>
/// The solution, either:
/// y = Solution[0]*x + Solution[1];
/// or if Solution.Length == 3:
/// y = Solution[0] + Solution[1]*x + Solution[2]*x*x;
/// </summary>
public double[] Solution { get; protected set; }
/// <summary>
/// Solve for y = mx + b ->
/// y = Solution[0]*x + Solution[1];
/// </summary>
/// <param name="x"></param>
/// <param name="y"></param>
public bool Solve_FirstOrder(double[] x, double[] y)
{
// This is the classical software approach to doing
// this fast, form the sums, do as little math as needed.
double sx = 0;
double sxx = 0;
double sy = 0;
double sxy = 0;
double syy = 0;
// compute the length.
int n = x.Length < y.Length ? x.Length : y.Length;
// Sum.
for (int i = 0; i < n; i++)
{
double xi = x[i];
double yi = y[i];
sx += xi;
sxx += xi * xi;
sy += yi;
sxy += xi * yi;
syy += yi * yi;
}
// y = bx + a
double b = ((n * sxy) - (sx * sy)) / ((n * sxx) - (sx * sx));
double a = (sy - (b*sx))/n;
// Store.
Solution = new double[2];
Solution[0] = a;
Solution[1] = b;
// Done.
return true;
}
/// <summary>
/// Solve for y = a0 + a1*x + a2*x^2 ->
/// y = Solution[0] + Solution[1]*x + Solution[2]*x*x;
/// </summary>
/// <param name="x"></param>
/// <param name="y"></param>
public bool Solve_SecondOrder(double[] x, double[] y)
{
// This is not done as a matrix operation
// because it is much faster to calculate this way in part
// because we know of the inherent matrix symetry and 3x3
// matrixes aren't soo bad to invert in code this way.
//
// Each is a sum:
// sx = sum of the x values.
// sx2 = sum of each x value after squaring that value.
double sx = 0;
double sx2 = 0;
double sx3 = 0;
double sx4 = 0;
double sy = 0;
double sxy = 0;
double sx2y = 0;
double sy2 = 0;
// n is the data set length.
int n = x.Length < y.Length ? x.Length : y.Length;
// Sum everything.
for (int i = 0; i < n; i++)
{
// Avoid repeated de-indexing, and repeated multiplies.
double xi = x[i];
double yi = y[i];
double xi2 = xi * xi;
// Sum.
sx += xi;
sx2 += xi2;
sy += yi;
sxy += xi * yi;
sy2 += yi * yi;
sx3 += xi2 * xi;
sx4 += xi2 * xi2;
sx2y += xi2 * yi;
}
// In matrix form, solving for [a0, a1, a2].
// n, sx, sx2 a0 sy
// sx, sx2, sx3 * a1 = sxy
// sx2, sx3, sx4 a2 sx2y
//
// Has a solution if the determinant isn't zero of:
// a, b, c sy a0
// b, d, e * sxy * 1/det = a1
// c, e, f sx2y a2
// det = n * a + sx * b + sx2 * c
double a = (sx2 * sx4) - (sx3 * sx3);
double b = (sx3 * sx2) - (sx * sx4);
double c = (sx * sx3) - (sx2 * sx2);
double d = (n * sx4) - (sx2 * sx2);
double e = (sx * sx2) - (n * sx3);
double f = (n * sx2) - (sx * sx);
double det = (n * a) + (sx * b) + (sx2 * c);
// Don't bother trying
if (det == 0) return false;
det = 1 / det;
double a0 = det * ((sy * a) + (sxy * b) + (sx2y * c));
double a1 = det * ((sy * b) + (sxy * d) + (sx2y * e));
double a2 = det * ((sy * c) + (sxy * e) + (sx2y * f));
// Store
Solution = new double[3];
Solution[0] = a0;
Solution[1] = a1;
Solution[2] = a2;
// Done.
return true;
}
/// <summary>
/// Compute the estimated Y values for this x.
/// </summary>
/// <param name="x"></param>
/// <returns></returns>
public double[] ComputeYValues(double[] x)
{
double[] y = new double[x.Length];
if (Solution.Length == 2)
{
double b = Solution[0];
double a = Solution[1];
for (int i = 0; i < x.Length; i++)
{
y[i] = (a * x[i]) + b;
}
return y;
}
if (Solution.Length == 3)
{
double a0 = Solution[0];
double a1 = Solution[1];
double a2 = Solution[2];
for (int i = 0; i < x.Length; i++)
{
double xi = x[i];
y[i] = a0 + (a1 * xi) + (a2 * xi * xi);
}
return y;
}
return null;
}
/// <summary>
/// Compute the error of this fit.
/// </summary>
/// <param name="x"></param>
/// <returns></returns>
public double ComputeError(double[] x, double[] y)
{
double[] ye = ComputeYValues(x);
double err = 0;
for (int i = 0; i < x.Length; i++)
{
double de = y[i] - ye[i];
err += (de * de);
}
err = System.Math.Sqrt(err) / x.Length;
return err;
}
}
}
|
By viewing downloads associated with this article you agree to the Terms of Service and the article's licence.
If a file you wish to view isn't highlighted, and is a text file (not binary), please
let us know and we'll add colourisation support for it.
Phil is a Principal Software developer focusing on weird yet practical algorithms that run the gamut of embedded and desktop (PID loops, Kalman filters, FFTs, client-server SOAP bindings, ASIC design, communication protocols, game engines, robotics).
In his personal life he is a part time mad scientist, full time dad, and studies small circle jujitsu, plays guitar and piano.