C#，数值计算——数据建模，直线拟合（Fitmed，最小化绝对偏差法）的计算方法与源程序

1 文本格式

using System;

namespace Legalsoft.Truffer
{
/// <summary>
/// fit a line minimizing absolute deviation
/// </summary>
public class Fitmed
{
private int ndata { get; set; }
private double a { get; set; }
private double b { get; set; }
private double abdev { get; set; }
private double[] x { get; set; }
private double[] y { get; set; }

public Fitmed(double[] xx, double[] yy)
{
this.ndata = xx.Length;
this.x = xx;
this.y = yy;

double sx = 0.0;
double sy = 0.0;
double sxy = 0.0;
double sxx = 0.0;
double chisq = 0.0;
for (int j = 0; j < ndata; j++)
{
sx += x[j];
sy += y[j];
sxy += x[j] * y[j];
sxx += Globals.SQR(x[j]);
}
double del = ndata * sxx - sx * sx;
a = (sxx * sy - sx * sxy) / del;
b = (ndata * sxy - sx * sy) / del;
for (int j = 0; j < ndata; j++)
{
double temp = y[j] - (a + b * x[j]);
chisq = temp * temp;
}
double sigb = Math.Sqrt(chisq / del);
double b1 = b;
double f1 = rofunc(b1);
if (sigb > 0.0)
{
double b2 = b + Globals.SIGN(3.0 * sigb, f1);
double f2 = rofunc(b2);
//if (b2 == b1)
if (Math.Abs(b2 - b1) <= float.Epsilon)
{
abdev /= ndata;
return;
}
while (f1 * f2 > 0.0)
{
b = b2 + 1.6 * (b2 - b1);
b1 = b2;
f1 = f2;
b2 = b;
f2 = rofunc(b2);
}
sigb = 0.01 * sigb;
while (Math.Abs(b2 - b1) > sigb)
{
b = b1 + 0.5 * (b2 - b1);
//if (b == b1 || b == b2)
if (Math.Abs(b - b1) <= float.Epsilon || Math.Abs(b - b2) <= float.Epsilon)
{
break;
}
double f = rofunc(b);
if (f * f1 >= 0.0)
{
f1 = f;
b1 = b;
}
else
{
f2 = f;
b2 = b;
}
}
}
abdev /= ndata;
}

public double rofunc(double b)
{
double EPS = float.Epsilon;
double[] arr = new double[ndata];
for (int j = 0; j < ndata; j++)
{
arr[j] = y[j] - b * x[j];
}
if ((ndata & 1) == 1)
{
a = Sorter.select((ndata - 1) >> 1, arr);
}
else
{
int jj = ndata >> 1;
a = 0.5 * (Sorter.select(jj - 1, arr) + Sorter.select(jj, arr));
}
abdev = 0.0;
double sum = 0.0;
for (int j = 0; j < ndata; j++)
{
double d = y[j] - (b * x[j] + a);
abdev += Math.Abs(d);
if (y[j] != 0.0)
{
d /= Math.Abs(y[j]);
}
if (Math.Abs(d) > EPS)
{
sum += (d >= 0.0 ? x[j] : -x[j]);
}
}
return sum;
}
}
}

2 代码格式

using System;namespace Legalsoft.Truffer
{/// <summary>/// fit a line minimizing absolute deviation/// </summary>public class Fitmed{private int ndata { get; set; }private double a { get; set; }private double b { get; set; }private double abdev { get; set; }private double[] x { get; set; }private double[] y { get; set; }public Fitmed(double[] xx, double[] yy){this.ndata = xx.Length;this.x = xx;this.y = yy;double sx = 0.0;double sy = 0.0;double sxy = 0.0;double sxx = 0.0;double chisq = 0.0;for (int j = 0; j < ndata; j++){sx += x[j];sy += y[j];sxy += x[j] * y[j];sxx += Globals.SQR(x[j]);}double del = ndata * sxx - sx * sx;a = (sxx * sy - sx * sxy) / del;b = (ndata * sxy - sx * sy) / del;for (int j = 0; j < ndata; j++){double temp = y[j] - (a + b * x[j]);chisq = temp * temp;}double sigb = Math.Sqrt(chisq / del);double b1 = b;double f1 = rofunc(b1);if (sigb > 0.0){double b2 = b + Globals.SIGN(3.0 * sigb, f1);double f2 = rofunc(b2);//if (b2 == b1)if (Math.Abs(b2 - b1) <= float.Epsilon){abdev /= ndata;return;}while (f1 * f2 > 0.0){b = b2 + 1.6 * (b2 - b1);b1 = b2;f1 = f2;b2 = b;f2 = rofunc(b2);}sigb = 0.01 * sigb;while (Math.Abs(b2 - b1) > sigb){b = b1 + 0.5 * (b2 - b1);//if (b == b1 || b == b2)if (Math.Abs(b - b1) <= float.Epsilon || Math.Abs(b - b2) <= float.Epsilon){break;}double f = rofunc(b);if (f * f1 >= 0.0){f1 = f;b1 = b;}else{f2 = f;b2 = b;}}}abdev /= ndata;}public double rofunc(double b){double EPS = float.Epsilon;double[] arr = new double[ndata];for (int j = 0; j < ndata; j++){arr[j] = y[j] - b * x[j];}if ((ndata & 1) == 1){a = Sorter.select((ndata - 1) >> 1, arr);}else{int jj = ndata >> 1;a = 0.5 * (Sorter.select(jj - 1, arr) + Sorter.select(jj, arr));}abdev = 0.0;double sum = 0.0;for (int j = 0; j < ndata; j++){double d = y[j] - (b * x[j] + a);abdev += Math.Abs(d);if (y[j] != 0.0){d /= Math.Abs(y[j]);}if (Math.Abs(d) > EPS){sum += (d >= 0.0 ? x[j] : -x[j]);}}return sum;}}
}