对于一般多项式:
K为多项式最高项次,a为不确定的常数项,共k+1个;
有离散数据集对应,其方差:
β为,方差函数S对β自变量第j个参数的梯度(偏导数):
当以上梯度为零时,S函数值最小,即:
中的每个每个偏导数构成一个等式:
...
则:
...
变为矩阵形式:
这样就变成线性方程求解形式,可用高斯消元等方法求得,注意在计算过程中要判断对角线上的值是否为零,如果等于零可以通过换行的方法解决;
/// <summary>
/// Function y = a0+a1*x+a2*x^2+ ... + an*x^n
/// </summary>
/// <param name="dataX">x values</param>
/// <param name="Y">y values</param>
/// <param name="N">Degree of a Polynomial</param>
/// <returns>[a0,a1,...,an]</returns>
public static double[] PolyRegress(double[] dataX, double[] Y,int N)
{
const double tiny = 0.00001;
int M = dataX.Length; // M = Number of data points
int D = N + ;
double x = , y = ;
double[,] V = new double[D, D];
double[] A = new double[D];
for (int i = ; i < M; i++)
{
x = dataX[i];
y = Y[i];
for(int j=;j<D;j++)
{
for(int k=;k<D;k++)
{
V[j, k] += Math.Pow(x, j + k);
}
A[j] += y * Math.Pow(x, j);
} }
for (int i = ; i < D; i++)
{
double m = V[i, i];
if (Math.Abs(m) < tiny)
{
for (int i2 = i + ; i2 < D; i2++)
{
if (Math.Abs(V[i2, i]) > tiny)
{
double tmp = ;
for (int c = ; c < D; c++)
{
tmp = V[i, c];
V[i, c] = V[i2, c];
V[i2, c] = tmp;
}
tmp = A[i];
A[i] = A[i2];
A[i2] = tmp;
break;
}
}
m = V[i, i];
} if (Math.Abs(m) > tiny)
{
for (int j = i; j < D; j++)
{
V[i, j] /= m;
}
A[i] /= m;
for (int k = ; k < D; k++)
{
if (k != i)
{
m = V[k, i];
for (int l = i; l < D; l++)
{
V[k, l] -= m * V[i, l];
}
A[k] -= m * A[i];
}
}
}
}
return A;
}