# LeastSquares.frink

``` use Matrix.frink /** This class finds least-squares curve fits.  It takes a set of points in x     and y and returns the best fit for a given curve type.     It allows you to specify the "basis functions" (or basis expressions) for     the curve type.  For example, if you wanted to find the best linear fit,     the basis functions would be  [x, 1].  For a quadratic (squared) fit, the     basis functions would be [x^2, x, 1].     This class finds the coefficients that best fit the provided points.  That     is, for the fit to a line mentioned above, this would calculate the     coefficients c1, c2 to best solve     y = c1 x + c2     This uses the Matrix.frink class to perform its solution, notably the     leastSquares[] method.     Curve-fitting can be performed on an overdetermined system, where there are     more measurements than equations.     This is the best discussion I've seen of least-squares fitting:     https://www.aleksandrhovhannisyan.com/blog/the-method-of-least-squares/     https://www.aleksandrhovhannisyan.com/blog/least-squares-fitting/     See:     https://mathworld.wolfram.com/LeastSquaresFitting.html     See also for special curve fits:     https://mathworld.wolfram.com/LeastSquaresFittingExponential.html     https://mathworld.wolfram.com/LeastSquaresFittingLogarithmic.html     https://mathworld.wolfram.com/LeastSquaresFittingPowerLaw.html */ class LeastSquares {    /** An array of x values to fit. */    var xvalues    /** An array of y values to fit. */    var yvalues    /** An array of basis expressions. */    var basisExprs    /** The variable which contains the curve fit coefficients.  It is a        one-column matrix */    var sol        new[xvalues, yvalues, basisExprs] :=    {       this.xvalues = xvalues       this.yvalues = yvalues       this.basisExprs = basisExprs       sol = fit[]    }    class fitLinear[xvalues, yvalues] :=    {       return fitDegree[xvalues, yvalues, 1]    }    class fitQuadratic[xvalues, yvalues] :=    {       return fitDegree[xvalues, yvalues, 2]    }    class fitCubic[xvalues, yvalues] :=    {       return fitDegree[xvalues, yvalues, 3]    }        class fitDegree[xvalues, yvalues, degree] :=    {       // Make basis functions like [x^2, x, 1]       a = new array       for i = 0 to degree           a.pushFirst[constructExpression["Power", [noEval[x], i]]]       return new LeastSquares[xvalues, yvalues, a]    }    /** Performs the internal curve fitting.  Solutions are placed into the        variable sol which is a one-column Matrix. */    fit[] :=    {       rows = length[xvalues]       cols = length[basisExprs]       a = new array[[rows,cols],0]       for row = 0 to rows-1          for col = 0 to cols-1          {             x = xvalues@row             a@row@col = eval[basisExprs@col]          }       A = new Matrix[a]       B = new Matrix[yvalues.transpose[]]       return A.leastSquares[B]    }    /** Returns an string representing the best fit.  For example, this        might return          "3.21 x + 1.54"        for a linear fit.    */    toExpressionString[] :=    {       cols = length[basisExprs]       estr = ""       for col = 0 to cols-1       {          estr = estr + "(" + inputForm[sol.get[col+1, 1]] + " * " + inputForm[basisExprs@col] + ")"          if col < cols-1             estr = estr + " + "       }              return estr    }    /** Returns an expression representing the best fit.  For example, this        might return          3.21 x + 1.54        for a linear fit.    */    toExpression[] :=    {       return parseToExpression[toExpressionString[]]    }    /** Returns an anonymous single-argument function that represents the best        fit, which you can then use to calculate additional y values. For        example, you could call it like:        f = toFunction[]        y1 = f        y2 = f    */    toFunction[] :=    {       return constructExpression["AnonymousFunction", [[x], toExpression[]]]    }        /** Returns the solution coefficients as a 1-column Matrix. */    toMatrix[] :=    {       return sol    }    /** Returns the solution coefficients as a row array. */    toArray[] :=    {       return sol.getColumnAsArray    }    /** Calculate the RMS of the residuals. */    residual[] :=    {       f = toFunction[]       size = length[xvalues]       r = 0 yvalues@0  // Make units work out       for i = 0 to size-1          r = r + (yvalues@i - f[xvalues@i])^2       return sqrt[r]    } } ```