RosettaCodeData/Task/Polynomial-regression/AWK/polynomial-regression.awk

85 lines
1.7 KiB
Awk

BEGIN{
i = 0;
xa[i] = 0; i++;
xa[i] = 1; i++;
xa[i] = 2; i++;
xa[i] = 3; i++;
xa[i] = 4; i++;
xa[i] = 5; i++;
xa[i] = 6; i++;
xa[i] = 7; i++;
xa[i] = 8; i++;
xa[i] = 9; i++;
xa[i] = 10; i++;
i = 0;
ya[i] = 1; i++;
ya[i] = 6; i++;
ya[i] = 17; i++;
ya[i] = 34; i++;
ya[i] = 57; i++;
ya[i] = 86; i++;
ya[i] =121; i++;
ya[i] =162; i++;
ya[i] =209; i++;
ya[i] =262; i++;
ya[i] =321; i++;
exit;
}
{
# (nothing to do)
}
END{
a = 0; b = 0; c = 0; # globals - will change by regression()
regression(xa,ya);
printf("y = %6.2f x^2 + %6.2f x + %6.2f\n",c,b,a);
printf("%-13s %-8s\n","Input","Approximation");
printf("%-6s %-6s %-8s\n","x","y","y^")
for (i=0;i<length(xa);i++) {
printf("%6.1f %6.1f %8.3f\n",xa[i],ya[i],eval(a,b,c,xa[i]));
}
}
function eval(a,b,c,x) {
return a+b*x+c*x*x;
}
# locals
function regression(x,y, n,xm,ym,x2m,x3m,x4m,xym,x2ym,sxx,sxy,sxx2,sx2x2,sx2y) {
n = 0
xm = 0.0;
ym = 0.0;
x2m = 0.0;
x3m = 0.0;
x4m = 0.0;
xym = 0.0;
x2ym = 0.0;
for (i in x) {
xm += x[i];
ym += y[i];
x2m += x[i] * x[i];
x3m += x[i] * x[i] * x[i];
x4m += x[i] * x[i] * x[i] * x[i];
xym += x[i] * y[i];
x2ym += x[i] * x[i] * y[i];
n++;
}
xm = xm / n;
ym = ym / n;
x2m = x2m / n;
x3m = x3m / n;
x4m = x4m / n;
xym = xym / n;
x2ym = x2ym / n;
sxx = x2m - xm * xm;
sxy = xym - xm * ym;
sxx2 = x3m - xm * x2m;
sx2x2 = x4m - x2m * x2m;
sx2y = x2ym - x2m * ym;
b = (sxy * sx2x2 - sx2y * sxx2) / (sxx * sx2x2 - sxx2 * sxx2);
c = (sx2y * sxx - sxy * sxx2) / (sxx * sx2x2 - sxx2 * sxx2);
a = ym - b * xm - c * x2m;
}