LeastSquares.frink

Download or view LeastSquares.frink in plain text format


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

   /** Construct and solve the system. */
   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[varname = "x"] :=
   {
      cols = length[basisExprs]
      estr = ""
      for col = 0 to cols-1
      {
         be2 = substituteExpression[basisExprs@col, parseToExpression["x"], parseToExpression[varname]]
         estr = estr + "(" + inputForm[sol.get[col+1, 1]] + " * " + inputForm[be2] + ")"
         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[varname = "x"] :=
   {
      return parseToExpression[toExpressionString[varname]]
   }

   /** 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[1]
       y2 = f[2]
   */

   toFunction[varname = "x"] :=
   {
      varsym = parseToExpression[varname]
      return constructExpression["AnonymousFunction", [[varsym], toExpression[varname]]]
   }

   
   /** Returns the solution coefficients as a 1-column Matrix. */
   toMatrix[] :=
   {
      return sol
   }

   /** Returns the solution coefficients as a row array. */
   toArray[] :=
   {
      return sol.getColumnAsArray[1]
   }

   /** Calculate the RMS of the residuals. */
   residual[] :=
   {
      f = toFunction[]
      size = length[xvalues]

      r = 0 (yvalues@0)^2  // Make units work out
      for i = 0 to size-1
         r = r + (yvalues@i - f[xvalues@i])^2

      return sqrt[r]
   }

   /** Calculates the r-value of the correlation between the yvalues and the fit
   function's predicted yvalues.  See
   https://stats.libretexts.org/Bookshelves/Introductory_Statistics/OpenIntro_Statistics_(Diez_et_al)./07%3A_Introduction_to_Linear_Regression/7.02%3A_Line_Fitting_Residuals_and_Correlation */

   rValue[] :=
   {
      f = toFunction[]
      ycalc = map[f, xvalues]
      len = length[yvalues]
      
      [meanx, sdx] = meanAndSD[yvalues, true]
      [meany, sdy] = meanAndSD[ycalc,   true]

      //println["mean x: $meanx  sdx: $sdx"]
      //println["mean y: $meany  sdy: $sdy"]
      sum = undef  // Make units come out right
      for i=rangeOf[yvalues]
      {
         term = ((yvalues@i - meanx) / sdx) * ((ycalc@i - meany) / sdy)
         if sum == undef
            sum = term
         else
            sum = sum + term
      }

      return sum / (len-1)
   }

   /** Calculates the mean and standard deviation of an array or
    enumerating expression.

    This uses Welford's algorithm as cited in Knuth, The Art of Computer
    Programming, Vol. 2, 3rd edition, page 232.

    Arguments:
     [list, sample]

    where
     list: is an array or enumerating expression of the items to average.
     sample: a boolean flag indicating if we want the standard deviation to be
         the sample standard variation (=true)  or the population standard
         variation (=false).

    Returns:
     [mean, sd, number]
   */

   class meanAndSD[list, sample] :=
   {
      M = 0 list@0     // Make units of measure come out right
      S = 0 (list@0)^2 // Make units of measure come out right
      k = 1

      for v = list
      {
         oldM = M
         diff = v-oldM
         M = M + diff / k
         S = S + diff * (v - M)
         k = k+1
      }

      if sample == true
         sub = 1            // Sample standard deviation, subtract 1 from num
      else
         sub = 0            // Population standard deviation
      
      return [M, sqrt[S/(k-1-sub)], k-1]
   }
}


Download or view LeastSquares.frink in plain text format


This is a program written in the programming language Frink.
For more information, view the Frink Documentation or see More Sample Frink Programs.

Alan Eliasen was born 19908 days, 4 hours, 39 minutes ago.