最小二乘回归二次曲线拟合
简介
最近的软件项目从一系列的数据点的二次曲线方程推导出的要求。这就是说，以确定A，B，C，Y = AX2 BX C.在确定A，B和C，我还需要R平方（决定系数）值。
快速搜索的谷歌未能带来一个合适的C＃类，很可能，这样做更彻底的搜索。
我有一个模糊的回忆称为"最小二乘回归"的东西，所以我去到Google。参考文献
最小二乘回归:{A}{A2}
R平方:{A3}
Cramer的规则:{A4}使用类
声明一个实例:LstSquQuadRegr solvr = new LstSquQuadRegr();
传递给它的一些点对（至少3个）:{C}
的值:double the_a_term = solvr.aTerm();
double the_b_term = solvr.bTerm();
double the_c_term = solvr.cTerm();
double the_rSquare = solvr.rSquare();
y = ax^2 + bx + c
，我们有一系列的点（X1，Y1），（X2，Y2）... ... （XN，YN）。for i = 1 to n
我们想要的值A，B，和C，最大限度地减少* XI ^ 2 bxi C.易偏差的平方的总和这种价值观的最佳拟合的二次方程。
让偏差的平方的总和是: n
F(a,b,c) = SUM (a*xi^2 + bxi + c - yi)^2.
i=1
dF/da = SUM 2*(a*xi^2+b*xi+c-yi)*xi^2 = 0,
dF/db = SUM 2*(a*xi^2+b*xi+c-yi)*xi = 0,
dF/dc = SUM 2*(a*xi^2+b*xi+c-yi) = 0.
（在这里，所有的款项范围内对i = 1，2，...，N）除以2并重新排列，使这三同时包含A，B和C三个未知数的线性方程组:(SUM xi^4)*a + (SUM xi^3)*b + (SUM xi^2)*c = SUM xi^2*yi,
(SUM xi^3)*a + (SUM xi^2)*b + (SUM xi)*c = SUM xi*yi,
(SUM xi^2)*a + (SUM xi)*b + (SUM 1)*c = SUM yi.
使用符号SJK意味着x_i ^ J * y_i ^ K的总和:a*S40 + b*S30 + c*S20 = S21
a*S30 + b*S20 + c*S10 = S11
a*S20 + b*S10 + c*S00 = S01
解联立方程组，利用Cramer的法律: [ S40 S30 S20 ] [ a ] [ S21 ]
[ S30 S20 S10 ] [ b ] = [ S11 ]
[ S20 S10 S00 ] [ c ] [ S01 ]
D = [ S40 S30 S20 ]
[ S30 S20 S10 ]
[ S20 S10 S00 ]
= S40(S20*S00 - S10*S10) - S30(S30*S00 - S10*S20) + S20(S30*S10 - S20*S20)
Da = [ S21 S30 S20 ]
[ S11 S20 S10 ]
[ S01 S10 S00 ]
= S21(S20*S00 - S10*S10) - S11(S30*S00 - S10*S20) + S01(S30*S10 - S20*S20)
Db = [ S40 S21 S20 ]
[ S30 S11 S10 ]
[ S20 S01 S00 ]
= S40(S11*S00 - S01*S10) - S30(S21*S00 - S01*S20) + S20(S21*S10 - S11*S20)
Dc = [ S40 S30 S21 ]
[ S30 S20 S11 ]
[ S20 S10 S01 ]
= S40(S20*S01 - S10*S11) - S30(S30*S01 - S10*S21) + S20(S30*S11 - S20*S21)
a = Da/D
b = Db/D
c = Dc/D
R2 = 1 - （剩余平方和/总平方和的总和）。 n
total sum of squares = SUM (yi - y_mean)^2.
i=1
这是测得的Y值和平均y值之间的差异的平方的总和。 n
residual sum of squares = SUM (yi - yi_predicted)^2.
i=1
这是测得的Y值和y的方程预测值之间的差异的平方的总和。"守则"
一堆辅助方法计算平方的各种款项。当计算A，B和C的值，我上面使用的，因为我发现它更容易跟踪SJK符号。/******************************************************************************
Class LstSquQuadRegr
A C# Class for Least Squares Regression for Quadratic Curve Fitting
Alex Etchells 2010
******************************************************************************/
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Collections;
public class LstSquQuadRegr
{
/* instance variables */
ArrayList pointArray = new ArrayList();
private int numOfEntries;
private double[] pointpair;
/*constructor */
public LstSquQuadRegr()
{
numOfEntries = 0;
pointpair = new double[2];
}
/*instance methods */
/// <summary>
/// add point pairs
/// </summary>
/// <param name="x">x value</param>
/// <param name="y">y value</param>
public void AddPoints(double x, double y)
{
pointpair = new double[2];
numOfEntries +=1;
pointpair[0] = x;
pointpair[1] = y;
pointArray.Add(pointpair);
}
/// <summary>
/// returns the a term of the equation ax^2 + bx + c
/// </summary>
/// <returns>a term</returns>
public double aTerm()
{
if (numOfEntries < 3)
{
throw new InvalidOperationException(
"Insufficient pairs of co-ordinates");
}
//notation sjk to mean the sum of x_i^j*y_i^k.
double s40 = getSx4(); //sum of x^4
double s30 = getSx3(); //sum of x^3
double s20 = getSx2(); //sum of x^2
double s10 = getSx(); //sum of x
double s00 = numOfEntries;
//sum of x^0 * y^0 ie 1 * number of entries
double s21 = getSx2y(); //sum of x^2*y
double s11 = getSxy(); //sum of x*y
double s01 = getSy(); //sum of y
//a = Da/D
return (s21*(s20 * s00 - s10 * s10) -
s11*(s30 * s00 - s10 * s20) +
s01*(s30 * s10 - s20 * s20))
/
(s40*(s20 * s00 - s10 * s10) -
s30*(s30 * s00 - s10 * s20) +
s20*(s30 * s10 - s20 * s20));
}
/// <summary>
/// returns the b term of the equation ax^2 + bx + c
/// </summary>
/// <returns>b term</returns>
public double bTerm()
{
if (numOfEntries < 3)
{
throw new InvalidOperationException(
"Insufficient pairs of co-ordinates");
}
//notation sjk to mean the sum of x_i^j*y_i^k.
double s40 = getSx4(); //sum of x^4
double s30 = getSx3(); //sum of x^3
double s20 = getSx2(); //sum of x^2
double s10 = getSx(); //sum of x
double s00 = numOfEntries;
//sum of x^0 * y^0 ie 1 * number of entries
double s21 = getSx2y(); //sum of x^2*y
double s11 = getSxy(); //sum of x*y
double s01 = getSy(); //sum of y
//b = Db/D
return (s40*(s11 * s00 - s01 * s10) -
s30*(s21 * s00 - s01 * s20) +
s20*(s21 * s10 - s11 * s20))
/
(s40 * (s20 * s00 - s10 * s10) -
s30 * (s30 * s00 - s10 * s20) +
s20 * (s30 * s10 - s20 * s20));
}
/// <summary>
/// returns the c term of the equation ax^2 + bx + c
/// </summary>
/// <returns>c term</returns>
public double cTerm()
{
if (numOfEntries < 3)
{
throw new InvalidOperationException(
"Insufficient pairs of co-ordinates");
}
//notation sjk to mean the sum of x_i^j*y_i^k.
double s40 = getSx4(); //sum of x^4
double s30 = getSx3(); //sum of x^3
double s20 = getSx2(); //sum of x^2
double s10 = getSx(); //sum of x
double s00 = numOfEntries;
//sum of x^0 * y^0 ie 1 * number of entries
double s21 = getSx2y(); //sum of x^2*y
double s11 = getSxy(); //sum of x*y
double s01 = getSy(); //sum of y
//c = Dc/D
return (s40*(s20 * s01 - s10 * s11) -
s30*(s30 * s01 - s10 * s21) +
s20*(s30 * s11 - s20 * s21))
/
(s40 * (s20 * s00 - s10 * s10) -
s30 * (s30 * s00 - s10 * s20) +
s20 * (s30 * s10 - s20 * s20));
}
public double rSquare() // get r-squared
{
if (numOfEntries < 3)
{
throw new InvalidOperationException(
"Insufficient pairs of co-ordinates");
}
// 1 - (residual sum of squares / total sum of squares)
return 1 - getSSerr() / getSStot();
}
/*helper methods*/
private double getSx() // get sum of x
{
double Sx = 0;
foreach (double[] ppair in pointArray)
{
Sx += ppair[0];
}
return Sx;
}
private double getSy() // get sum of y
{
double Sy = 0;
foreach (double[] ppair in pointArray)
{
Sy += ppair[1];
}
return Sy;
}
private double getSx2() // get sum of x^2
{
double Sx2 = 0;
foreach (double[] ppair in pointArray)
{
Sx2 += Math.Pow(ppair[0], 2); // sum of x^2
}
return Sx2;
}
private double getSx3() // get sum of x^3
{
double Sx3 = 0;
foreach (double[] ppair in pointArray)
{
Sx3 += Math.Pow(ppair[0], 3); // sum of x^3
}
return Sx3;
}
private double getSx4() // get sum of x^4
{
double Sx4 = 0;
foreach (double[] ppair in pointArray)
{
Sx4 += Math.Pow(ppair[0], 4); // sum of x^4
}
return Sx4;
}
private double getSxy() // get sum of x*y
{
double Sxy = 0;
foreach (double[] ppair in pointArray)
{
Sxy += ppair[0] * ppair[1]; // sum of x*y
}
return Sxy;
}
private double getSx2y() // get sum of x^2*y
{
double Sx2y = 0;
foreach (double[] ppair in pointArray)
{
Sx2y += Math.Pow(ppair[0], 2) * ppair[1]; // sum of x^2*y
}
return Sx2y;
}
private double getYMean() // mean value of y
{
double y_tot = 0;
foreach (double[] ppair in pointArray)
{
y_tot += ppair[1];
}
return y_tot/numOfEntries;
}
private double getSStot() // total sum of squares
{
//the sum of the squares of the differences between
//the measured y values and the mean y value
double ss_tot = 0;
foreach (double[] ppair in pointArray)
{
ss_tot += Math.Pow(ppair[1] - getYMean(), 2);
}
return ss_tot;
}
private double getSSerr() // residual sum of squares
{
//the sum of the squares of te difference between
//the measured y values and the values of y predicted by the equation
double ss_err = 0;
foreach (double[] ppair in pointArray)
{
ss_err += Math.Pow(ppair[1] - getPredictedY(ppair[0]), 2);
}
return ss_err;
}
private double getPredictedY(double x)
{
//returns value of y predicted by the equation for a given value of x
return aTerm() * Math.Pow(x, 2) + bTerm() * x + cTerm();
}
}
这是真的 - 它似乎同意由Excel给出的值非常接近。我希望这样可以节省别人的麻烦工作。
非常干净的解决方案。工程完全在我的应用程序。感谢这么多救了我重新发明轮子的麻烦。它的行为精美的抖动过滤器