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 as
[x, y] pairs and returns the best fit for a given curve type.
Note that if you have x and y points in separate arrays, you should call
zip[x,y] to turn them into an array of [x,y] pairs before passing in data.
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.
See the leastSquaresTest.frink program for examples of use of this library.
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
Also see curveFit.frink for clear examples of linear and quadratic fit.
*/
class LeastSquares
{
/** An array of [x,y] pairs to fit. */
var data
/** 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. data is an array or enumerating
expression of [x,y] pairs. */
new[data, basisExprs] :=
{
this.data = toArray[data]
this.basisExprs = basisExprs
sol = fit[]
}
/** Class-level method to fit a dataset to a line.
data is an array or enumerating expression of [x,y] pairs. */
class fitLinear[data] :=
{
return fitDegree[data, 1]
}
/** Class-level method to fit a dataset to a quadratic (squared) function.
data is an array or enumerating expression of [x,y] pairs. */
class fitQuadratic[data] :=
{
return fitDegree[data, 2]
}
/** Class-level method to fit a dataset to a cubic function.
data is an array or enumerating expression of [x,y] pairs. */
class fitCubic[data] :=
{
return fitDegree[data, 3]
}
/** Class-level method to fit a dataset to a specified degree polynomial.
data is an array or enumerating expression of [x,y] pairs. */
class fitDegree[data, degree] :=
{
// Make basis functions like [x^2, x, 1]
a = new array
for i = 0 to degree
a.pushFirst[constructExpression["Power", [noEval[x], i]]]
// The noEval prevents the constructor from trying to eval and throwing
// a bunch of "unknown variable" messages.
return new LeastSquares[data, noEval[a]]
}
/** Fits a polynomial of the specified degree in 2 pieces. There should be
at least 2*(degree+1) elements in the input list.
data is an array or enumerating expression of [x,y] pairs.
returns:
[func, rvalue]
*/
class fitDegreePiecewise[data, degree] :=
{
sorted = sort[toArray[data], byColumn[0]]
len = length[sorted]
bestR = -1000
bestF = undef
for cutBefore = degree+1 to len-(degree+1)
{
cutx = sorted@cutBefore@0 // x value after cut point
data1 = first[sorted, cutBefore]
data2 = rest[sorted, cutBefore]
e1 = fitDegree[data1, degree].toExpression[]
e2 = fitDegree[data2, degree].toExpression[]
fs = "{|x| \n" +
" if x < " + inputForm[cutx] + "\n" +
" return " + inputForm[e1] + "\n" +
" else\n" +
" return " + inputForm[e2] + "\n" +
"}"
if debug
println[fs]
f = parseToExpression[fs]
r = rValue[f, sorted]
if debug
print["r = $r "]
if r > bestR
{
bestR = r
bestF = f
if debug
println["*"]
}
if debug
println["\n"]
}
return [bestF, bestR]
}
/** Performs the internal curve fitting. Solutions are placed into the
variable sol which is a one-column Matrix. */
fit[] :=
{
rows = length[data]
cols = length[basisExprs]
a = new array[[rows,cols],0]
for row = 0 to rows-1
for col = 0 to cols-1
{
[x,y] = data@row
a@row@col = eval[basisExprs@col]
}
A = new Matrix[a]
B = new Matrix[data.getColumn[1].transpose[]] // Tranpose y values
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[data]
r = 0 (data@0@1)^2 // Make units work out (this is first y value)
for i = 0 to size-1
{
[x,y] = data@i
r = r + (y - f[x])^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[]
xvalues = data.getColumn[0]
yvalues = data.getColumn[1]
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)
}
}
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 20217 days, 15 hours, 25 minutes ago.