/** This library contains routines to perform a best linear or other type of fit to a set of data points. This is often referred to as "regression" because of some historical accident. These functions are designed to preserve units of measure, to provide coefficients with correct units of measure, and to work for purely symbolic values! To use this file, see the examples in curveFitTest.frink Those examples are powerful and let you, say, derive the gravitic equation. Also see LeastSquares.frink and the test program LeastSquaresTest.frink which allow you to fit arbitrary linear systems to arbitrary basis functions. Also see the leastSquares function in Matrix.frink and the test program MatrixQRTest.frink which demonstrate linear least-squares fitting using QR decomposition which is general for a wide variety of linear systems with non-exact measurements. */ /** Performs a best linear fit of the specified data points. In other words, this finds the coefficients a and b of a line (in slope-intercept form) with the equation: y = a x + b params: data: an array or set of [x,y] pairs. returns: [a, b, r] where r is the correlation coefficient */ linearFit[data] := { array = toArray[data] N = length[array] if N < 2 return [undef, undef, undef] // We do it this way to preserve units of measure by initializing the sum // with the first element. [x,y] = array@0 sxy = x y // Sum of x*y sx = x // Sum of x sy = y // Sum of y sx2 = x^2 // Sum of x^2 sy2 = y^2 // Sum of y^2 for i = 1 to N-1 { [x,y] = array@i sxy = sxy + x y sx = sx + x sy = sy + y sx2 = sx2 + x^2 sy2 = sy2 + y^2 } denom = N sx2 - sx^2 a = (N sxy - sx sy) / denom b = (sy sx2 - sx sxy) / denom // println["got here, denom=$denom"] // println["N=$N, sxy=$sxy, sx=$sx, sy=$sy"] // TODO: make these arbitrary-precision sqrt? That means we have to // fix arbitrary-precision sqrt so it takes units of measure. r = (N sxy - sx sy) / (sqrt[denom] sqrt[N sy2 - sy^2]) // println["Got past r"] return [a, b, r] } /** This performs a linear fit of the specified data as above, but instead of returning the coefficients, it returns a function representing the line, or, in other words, an anonymous function f that represents y = f[x] of the best-fit line, or in more Frink-like notation, the anonymous function looks like: {|x| a x + b } */ linearFitFunction[data] := { [a, b, r] = linearFit[data] return eval["{|x| (" + inputForm[a] + ") x + (" + inputForm[b] + ")}"] } /** Performs a best quadratic fit of the specified data points. In other words, this finds the coefficients a,b,c of a quadratic equation, that is, the equation: y = a x^2 + b x + c params: data: an array or set of [x,y] pairs. returns: [a, b, c, r] where r is the correlation coefficient */ quadraticFit[data] := { array = toArray[data] N = length[array] if N < 3 return [undef, undef, undef, undef] // We do it this way to preserve units of measure by initializing the sum // with the first element. [x,y] = array@0 sx = x // Sum of x sy = y // Sum of y sxy = x y // Sum of x*y sx2 = x^2 // Sum of x^2 sx3 = x^3 // Sum of x^3 sx4 = x^4 // Sum of x^4 sy2 = y^2 // Sum of y^2 sx2y = x^2 y // Sum of x^2 * y for i = 1 to N-1 { [x,y] = array@i sxy = sxy + x y sx = sx + x sy = sy + y sx2 = sx2 + x^2 sx3 = sx3 + x^3 sx4 = sx4 + x^4 sy2 = sy2 + y^2 sx2y = sx2y + x^2 y } // These symbol changes make it more concise and match the symbols // in Jean Meeus, Astronomical Algorithms, 4.5 and 4.6 P = sx Q = sx2 R = sx3 S = sx4 T = sy U = sxy V = sx2y D = N Q S + 2 P Q R - Q^3 - P^2 S - N R^2 a = (N Q V + P R T + P Q U - Q^2 T - P^2 V - N R U) / D b = (N S U + P Q V + Q R T - Q^2 U - P S T - N R V) / D c = (Q S T + Q R U + P R V - Q^2 V - P S U - R^2 T) / D meany = sy / N [x,y] = array@0 SSE = (y - a x^2 - b x - c)^2 SST = (y - meany)^2 for i = 1 to N-1 { [x,y] = array@i SSE = SSE + (y - a x^2 - b x - c)^2 SST = SST + (y - meany)^2 } r = sqrt[1 - SSE/SST] return [a, b, c, r] } /** This performs a quadratic fit of the specified data as above, but instead of returning the coefficients, it returns a function representing the function, or, in other words, an anonymous function f that represents y = f[x] of the best-fit line, or in more Frink-like notation, the anonymous function looks like: {|x| a x^2 + b x + c } */ quadraticFitFunction[data] := { [a, b, c, r] = quadraticFit[data] return eval["{|x| (" + inputForm[a] + ") x^2 + (" + inputForm[b] + ") x + (" + inputForm[c] + ")}"] }