C# 最小二乘法拟合
最小二乘法（又称最小平方法）是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。
利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。
public enum RegressionLine
{
kNone = 0,
kLiner = 1, // 线性
kPower = 2, // 幂次
kLogarithm = 3, // 对数
kExponent = 4, // 指数
kMultiLinear = 5, // 多项式
}
public abstract class RegressionAlgo
{
public virtual RegressionLine GetCategory() { return RegressionLine.kNone; }
protected double[] OriginX = null;
protected double[] OriginY = null;
// 设置点集
// 外面应该剔除无效值以后再传进来,里面不做判断
public virtual void SetOriginData(double[] x, double[] y)
{
this.OriginX = x;
this.OriginY = y;
}
// 计算回归线参数
public abstract bool CalculateRegression();
// 根据X计算Y
public abstract double CalculateValue(double x);
// 获取参数列表
public abstract double[] GetResult();
// 获取R2
public abstract double GetR2();
// 获取公式字符串
public abstract string GetFormula();
}
public class ShRegLinear : RegressionAlgo
{
public override RegressionLine GetCategory() { return RegressionLine.kLiner; }
public double A
{
get;
private set;
}
public double B
{
get;
private set;
}
public double R2
{
get;
private set;
}
public ShRegLinear()
{
// Y = a * X + b
this.A = 0.0;
this.B = 0.0;
}
public override string GetFormula()
{
return string.Format("Y = {0:F4} * X + {1:F4} \r\n R2 ={2:F4}", this.A, this.B, this.R2);
}
public ShRegLinear(double[] param, double r2)
{
if (param.Length == 2)
{
this.A = param[0];
this.B = param[1];
}
this.R2 = r2;
}
public override bool CalculateRegression()
{
if (this.OriginX == null || this.OriginX.Length < 3 ||
this.OriginY == null || this.OriginX.Length != this.OriginY.Length)
return false;
// 最小二乘法 https://zh.wikipedia.org/wiki/%E6%9C%80%E5%B0%8F%E4%BA%8C%E4%B9%98%E6%B3%95
double averX = this.OriginX.Average();
double averY = this.OriginY.Average();
double differencexy = 0.0;
double sumxy = 0.0;
double sumx = 0.0;
double sumy = 0.0;
for (int i = 0; i < this.OriginX.Length; i++)
{
differencexy += (this.OriginX[i] - averX) * (this.OriginY[i] - averY);
sumxy += this.OriginX[i] * this.OriginY[i];
sumx += Math.Pow((this.OriginX[i] - averX), 2);
sumy += Math.Pow((this.OriginY[i] - averY), 2);
}
// 经验回归系数
this.A = differencexy / sumx;
this.B = averY - this.A * averX;
this.R2 = Math.Pow(differencexy, 2) / (sumx * sumy);
return true;
}
public override double CalculateValue(double x)
{
return (this.A * x + this.B);
}
public override double[] GetResult()
{
double[] result = new double[]
{
this.A,
this.B,
};
return result;
}
public override double GetR2()
{
return this.R2;
}
}
public class ShRegExponential : RegressionAlgo
{
public override RegressionLine GetCategory() { return RegressionLine.kExponent; }
public double A
{
get
{
return Math.Exp(this.RegLinear.A);
}
}
public double B
{
get
{
return Math.Exp(this.RegLinear.B);
}
}
public double R2
{
get
{
return this.RegLinear.R2;
}
}
private ShRegLinear RegLinear = null;
public ShRegExponential()
{
// Y = b * a ^ X
this.RegLinear = new ShRegLinear();
}
public override string GetFormula()
{
return string.Format("Y = {1:F4} * {0:F4} ^ X \r\nR2 = {2:F4}", this.A, this.B, this.R2);
}
public ShRegExponential(double[] param, double r2)
{
if (param.Length == 2)
{
double[] temp = new double[2];
temp[0] = Math.Log(param[0]);
temp[1] = Math.Log(param[1]);
this.RegLinear = new ShRegLinear(temp, r2);
}
}
public override bool CalculateRegression()
{
if (this.OriginX == null || this.OriginX.Length < 3 ||
this.OriginY == null || this.OriginX.Length != this.OriginY.Length)
return false;
List<double> x = new List<double>();
List<double> y = new List<double>();
for (int i = 0; i < this.OriginY.Length; i++)
{
if (this.OriginY[i] > 0)
{
x.Add(this.OriginX[i]);
y.Add(Math.Log(this.OriginY[i]));
}
}
this.RegLinear.SetOriginData(x.ToArray(), y.ToArray());
return this.RegLinear.CalculateRegression();
}
public override double CalculateValue(double x)
{
return (this.B * Math.Pow(this.A, x));
}
public override double[] GetResult()
{
double[] result = new double[]
{
this.A,
this.B,
};
return result;
}
public override double GetR2()
{
return this.R2;
}
}
public class ShRegPower : RegressionAlgo
{
public override RegressionLine GetCategory() { return RegressionLine.kPower; }
public double A
{
get
{
return this.RegLinear.A;
}
}
public double B
{
get
{
return Math.Exp(this.RegLinear.B);
}
}
public double R2
{
get
{
return this.RegLinear.R2;
}
}
private ShRegLinear RegLinear = null;
public ShRegPower()
{
// Y = x^a * b
this.RegLinear = new ShRegLinear();
}
public override string GetFormula()
{
return string.Format("Y = X ^ {0:F4} * {1:F4} \r\nR2 = {2:F4}", this.A, this.B, this.R2);
}
public ShRegPower(double[] param, double r2)
{
if (param.Length == 2)
{
double[] temp = new double[2];
temp[0] = param[0];
temp[1] = Math.Log(param[1]);
this.RegLinear = new ShRegLinear(temp, r2);
}
}
public override bool CalculateRegression()
{
if (this.OriginX == null || this.OriginX.Length < 3 ||
this.OriginY == null || this.OriginX.Length != this.OriginY.Length)
return false;
List<double> x = new List<double>();
List<double> y = new List<double>();
for (int i = 0; i < this.OriginX.Length; i++)
{
if (this.OriginX[i] > 0 && this.OriginY[i] > 0)
{
x.Add(Math.Log(this.OriginX[i]));
y.Add(Math.Log(this.OriginY[i]));
}
}
this.RegLinear.SetOriginData(x.ToArray(), y.ToArray());
return this.RegLinear.CalculateRegression();
}
public override double CalculateValue(double x)
{
return this.B * Math.Pow(x, this.A);
}
public override double[] GetResult()
{
double[] result = new double[]
{
this.A,
this.B,
};
return result;
}
public override double GetR2()
{
return this.R2;
}
}
public class ShRegLogarithm : RegressionAlgo
{
public override RegressionLine GetCategory() { return RegressionLine.kLogarithm; }
public double A
{
get
{
return this.RegLinear.A;
}
}
public double B
{
get
{
return this.RegLinear.B;
}
}
public double R2
{
get
{
return this.RegLinear.R2;
}
}
private ShRegLinear RegLinear = null;
public ShRegLogarithm()
{
// Y = a *(ln x) + b
this.RegLinear = new ShRegLinear();
}
public override string GetFormula()
{
return string.Format("Y = {0:F4} * ln X + {1:F4} \r\nR2 = {2:F4}", this.A, this.B, this.R2);
}
public ShRegLogarithm(double[] param, double r2)
{
this.RegLinear = new ShRegLinear(param, r2);
}
public override bool CalculateRegression()
{
if (this.OriginX == null || this.OriginX.Length < 3 ||
this.OriginY == null || this.OriginX.Length != this.OriginY.Length)
return false;
List<double> x = new List<double>();
List<double> y = new List<double>();
for (int i = 0; i < this.OriginX.Length; i++)
{
if (this.OriginX[i] > 0)
{
x.Add(Math.Log(this.OriginX[i]));
y.Add(this.OriginY[i]);
}
}
this.RegLinear.SetOriginData(x.ToArray(), y.ToArray());
return this.RegLinear.CalculateRegression();
}
public override double CalculateValue(double x)
{
return this.A * Math.Log(x) + this.B;
}
public override double[] GetResult()
{
double[] result = new double[]
{
this.A,
this.B,
};
return result;
}
public override double GetR2()
{
return this.R2;
}
}
public class ShRegMultiLinear : RegressionAlgo
{
public override RegressionLine GetCategory() { return RegressionLine.kMultiLinear; }
public double[] Param = null;
public double R2
{
get;
private set;
}
public ShRegMultiLinear()
{
// Y = a[0] + a[1] * X + a[2] * X + ......
}
public override string GetFormula()
{
string text = "Y = ";
for (int index = this.Param.Length - 1; index >= 0; index--)
{
if (Math.Abs(this.Param[index]) > 1e-5)
{
if (index > 0)
{
text += string.Format("{0} * X^{1}", this.Param[index], index);
}
else
{
text += string.Format("{0}", this.Param[index]);
}
}
}
text += string.Format("\r\nR2 = {0}", this.R2);
return text;
}
public ShRegMultiLinear(double[] param, double r2)
{
this.Param = param;
this.R2 = r2;
}
public void SetDegree(int ex)
{
// 多项式回归至少是2次
if (ex < 2)
ex = 2;
this.Param = new double[ex + 1];
}
public override bool CalculateRegression()
{
// 几次多项式则至少要有这么多点
if (this.OriginX == null || this.Param == null || this.OriginX.Length < this.Param.Length ||
this.OriginY == null || this.OriginX.Length != this.OriginY.Length)
return false;
// 计算方程组的增广矩阵
double[,] em = new double[this.Param.Length, this.Param.Length + 1];
for (int i = 0; i < this.Param.Length; i++)
{
for (int j = 0; j < this.Param.Length; j++)
{
for (int k = 0; k < this.OriginX.Length; k++)
{
em[i, j] += Math.Pow(this.OriginX[k], i + j);
}
}
for (int k = 0; k < this.OriginX.Length; k++)
{
em[i, this.Param.Length] += this.OriginY[k] * Math.Pow(this.OriginX[k], i);
}
}
// 求解方程(消元过程)
for (int k = 0; k < this.Param.Length - 1; k++)
{
for (int i = k + 1; i < this.Param.Length; i++)
{
double t = -em[i, k] / em[k, k];
for (int j = k + 1; j <= this.Param.Length; j++)
{
em[i, j] += em[k, j] * t;
}
}
}
// 回代求解
for (int i = this.Param.Length - 1; i >= 0; i--)
{
for (int j = i + 1; j < this.Param.Length; j++)
{
em[i, this.Param.Length] -= em[i, j] * em[j, this.Param.Length];
}
em[i, this.Param.Length] /= em[i, i];
}
// 获取参数
this.Param[0] = em[0, this.Param.Length];
for (int i = 1; i < this.Param.Length; i++)
{
this.Param[i] = em[i, this.Param.Length];
}
return true;
}
public override double CalculateValue(double x)
{
if (this.Param.Length < 3)
return DefaultField.InvalidValue;
double sum = 0.0;
for (int index = 0; index < this.Param.Length; index++)
{
sum += this.Param[index] * Math.Pow(x, index);
}
return sum;
}
public override double[] GetResult()
{
return this.Param;
}
public override double GetR2()
{
double averY = this.OriginY.Average();
double ssr = 0.0;
double sse = 0.0;
for (int i = 0; i < this.OriginX.Length; i++)
{
double yi = CalculateValue(this.OriginX[i]);
// 回归平方和
ssr += Math.Pow(yi - averY, 2);
// 残差平方和
sse += Math.Pow(yi - this.OriginY[i], 2);
}
return 1 - (sse / (ssr + sse));
}
}