读书人

最小二乘法算法有关问题

发布时间: 2012-12-27 10:17:10 作者: rapoo

最小二乘法算法问题。
我做了一个最小二乘法的算法,是用C#写的。但是这个算法是有截距的,我希望截距等于0.求高人指教如何修改,或者提供一个新的算法,不要截距。类似y=ax2+bx。


///<summary>
///用最小二乘法拟合二元多次曲线
///例如y=ax+b
///其中MultiLine将返回a,b两个参数。
///a对应MultiLine[1]
///b对应MultiLine[0]
///</summary>
///<param name="arrX">已知点的x坐标集合</param>
///<param name="arrY">已知点的y坐标集合</param>
///<param name="length">已知点的个数</param>
///<param name="dimension">方程的最高次数</param>
public static double[] MultiLine(double[] arrX, double[] arrY, int length, int dimension)//二元多次线性方程拟合曲线
{
int n = dimension+1; //dimension次方程需要求 dimension+1个 系数
double[,] Guass = new double[n, n+1]; //高斯矩阵 例如:y=a0+a1*x+a2*x*x
for (int i = 0; i < n; i++)
{
int j;
for (j = 0; j < n; j++)
{
Guass[i, j] = SumArr(arrX, j + i, length);
}
Guass[i, j] = SumArr(arrX, i, arrY, 1, length);
}

return ComputGauss(Guass, n);

}
private static double SumArr(double[] arr, int n, int length) //求数组的元素的n次方的和
{
double s = 0;
for (int i = 0; i < length; i++)
{
if (arr[i] != 0 || n != 0)
s = s + Math.Pow(arr[i], n);
else


s = s + 1;
}
return s;
}
private static double SumArr(double[] arr1, int n1, double[] arr2, int n2, int length)
{
double s = 0;
for (int i = 0; i < length; i++)
{
if ((arr1[i] != 0 || n1 != 0) && (arr2[i] != 0 || n2 != 0))
s = s + Math.Pow(arr1[i], n1) * Math.Pow(arr2[i], n2);
else
s = s + 1;
}
return s;

}
private static double[] ComputGauss(double[,] Guass, int n)
{
int i, j;
int k, m;
double temp;
double max;
double s;
double[] x = new double[n];

for (i = 0; i < n; i++) x[i] = 0.0;//初始化


for (j = 0; j < n; j++)
{
max = 0;

k = j;
for (i = j; i < n; i++)
{
if (Math.Abs(Guass[i, j]) > max)
{
max = Guass[i, j];


k = i;
}
}



if (k != j)
{
for (m = j; m < n + 1; m++)
{
temp = Guass[j, m];
Guass[j, m] = Guass[k, m];
Guass[k, m] = temp;

}
}

if (0 == max)
{
// "此线性方程为奇异线性方程"

return x;
}


for (i = j + 1; i < n; i++)
{

s = Guass[i, j];
for (m = j; m < n + 1; m++)
{
Guass[i, m] = Guass[i, m] - Guass[j, m] * s / (Guass[j, j]);

}
}


}//结束for (j=0;j<n;j++)


for (i = n - 1; i >=0; i--)
{


s = 0;
for (j = i + 1; j < n; j++)
{
s = s + Guass[i, j] * x[j];
}

x[i] = (Guass[i, n] - s) / Guass[i, i];

}

return x;
}
}
}//返回值是函数的系数


[解决办法]
只要简单修改就可以了
思路:让最低次的系数,不参与循环,直接使用默认值0.
改两个地方:
1. MultiLine第六行: for (j = 0; j < n; j++) 改为 for (j = 1; j < n; j++)
2. ComputGauss第八行(不计算空行): for (j = 0; j < n; j++) 改为 for (j = 1; j < n; j++)

读书人网 >软件架构设计

热点推荐