Polynomial regression: Difference between revisions
(→{{header|Python}}: ++ gnu octave) |
|||
Line 52: | Line 52: | ||
</pre> |
</pre> |
||
=={{header|ALGOL 68}}== |
|||
{{trans|Ada}} |
|||
{{works with|ALGOL 68|Standard - ''lu decomp'' and ''lu solve'' are from the [[GSL]] library}} |
|||
{{works with|ALGOL 68G|Any - tested with release mk15-0.8b.fc9.i386}} |
|||
<!-- {{does not work with|ELLA ALGOL 68|Any (with appropriate job cards AND formatted transput statements removed) - tested with release 1.8.8d.fc9.i386 - ELLA has no FORMATted transput}} --> |
|||
<lang algol>MODE FIELD = REAL; |
|||
MODE |
|||
VEC = [0]FIELD, |
|||
MAT = [0,0]FIELD; |
|||
PROC VOID raise index error := VOID: ( |
|||
print(("stop", new line)); |
|||
stop |
|||
); |
|||
COMMENT from http://rosettacode.org/wiki/Matrix_Transpose#ALGOL_68 END COMMENT |
|||
OP ZIP = ([,]FIELD in)[,]FIELD:( |
|||
[2 LWB in:2 UPB in,1 LWB in:1UPB in]FIELD out; |
|||
FOR i FROM LWB in TO UPB in DO |
|||
out[,i]:=in[i,] |
|||
OD; |
|||
out |
|||
); |
|||
COMMENT http://rosettacode.org/wiki/Matrix_multiplication#ALGOL_68 END COMMENT |
|||
OP * = (VEC a,b)FIELD: ( # basically the dot product # |
|||
FIELD result:=0; |
|||
IF LWB a/=LWB b OR UPB a/=UPB b THEN raise index error FI; |
|||
FOR i FROM LWB a TO UPB a DO result+:= a[i]*b[i] OD; |
|||
result |
|||
); |
|||
OP * = (VEC a, MAT b)VEC: ( # overload vector times matrix # |
|||
[2 LWB b:2 UPB b]FIELD result; |
|||
IF LWB a/=LWB b OR UPB a/=UPB b THEN raise index error FI; |
|||
FOR j FROM 2 LWB b TO 2 UPB b DO result[j]:=a*b[,j] OD; |
|||
result |
|||
); |
|||
OP * = (MAT a, b)MAT: ( # overload matrix times matrix # |
|||
[LWB a:UPB a, 2 LWB b:2 UPB b]FIELD result; |
|||
IF 2 LWB a/=LWB b OR 2 UPB a/=UPB b THEN raise index error FI; |
|||
FOR k FROM LWB result TO UPB result DO result[k,]:=a[k,]*b OD; |
|||
result |
|||
); |
|||
COMMENT from http://rosettacode.org/wiki/Pyramid_of_numbers#ALGOL_68 END COMMENT |
|||
OP / = (VEC a, MAT b)VEC: ( # vector division # |
|||
[LWB a:UPB a,1]FIELD transpose a; |
|||
transpose a[,1]:=a; |
|||
(transpose a/b)[,1] |
|||
); |
|||
OP / = (MAT a, MAT b)MAT:( # matrix division # |
|||
[LWB b:UPB b]INT p ; |
|||
INT sign; |
|||
[,]FIELD lu = lu decomp(b, p, sign); |
|||
[LWB a:UPB a, 2 LWB a:2 UPB a]FIELD out; |
|||
FOR col FROM 2 LWB a TO 2 UPB a DO |
|||
out[,col] := lu solve(b, lu, p, a[,col]) [@LWB out[,col]] |
|||
OD; |
|||
out |
|||
); |
|||
FORMAT int repr = $g(0)$, |
|||
real repr = $g(-7,4)$; |
|||
PROC fit = (VEC x, y, INT order)VEC: |
|||
BEGIN |
|||
[0:order, LWB x:UPB x]FIELD a; # the plane # |
|||
FOR i FROM 2 LWB a TO 2 UPB a DO |
|||
FOR j FROM LWB a TO UPB a DO |
|||
a [j, i] := x [i]**j |
|||
OD |
|||
OD; |
|||
( y * ZIP a ) / ( a * ZIP a ) |
|||
END # fit #; |
|||
PROC print polynomial = (VEC x)VOID: ( |
|||
BOOL empty := TRUE; |
|||
FOR i FROM UPB x BY -1 TO LWB x DO |
|||
IF x[i] NE 0 THEN |
|||
IF x[i] > 0 AND NOT empty THEN print ("+") FI; |
|||
empty := FALSE; |
|||
IF x[i] NE 1 OR i=0 THEN |
|||
IF ENTIER x[i] = x[i] THEN |
|||
printf((int repr, x[i])) |
|||
ELSE |
|||
printf((real repr, x[i])) |
|||
FI |
|||
FI; |
|||
CASE i+1 IN |
|||
SKIP,print(("x")) |
|||
OUT |
|||
printf(($"x**"g(0)$,i)) |
|||
ESAC |
|||
FI |
|||
OD; |
|||
IF empty THEN print("0") FI; |
|||
print(new line) |
|||
); |
|||
fitting: BEGIN |
|||
VEC c = |
|||
fit |
|||
( (0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0), |
|||
(1.0, 6.0, 17.0, 34.0, 57.0, 86.0, 121.0, 162.0, 209.0, 262.0, 321.0), |
|||
2 |
|||
); |
|||
print polynomial(c); |
|||
VEC d = |
|||
fit |
|||
( (0, 1, 2, 3, 4, 5, 6, 7, 8, 9), |
|||
(2.7, 2.8, 31.4, 38.1, 58.0, 76.2, 100.5, 130.0, 149.3, 180.0), |
|||
2 |
|||
); |
|||
print polynomial(d) |
|||
END # fitting #</lang> |
|||
output: |
|||
<pre> |
|||
3x**2+2x+1 |
|||
1.0848x**2+10.3552x-0.6164 |
|||
</pre> |
|||
=={{header|C}}== |
=={{header|C}}== |
||
{{libheader|libgsl}} |
{{libheader|libgsl}} |
Revision as of 13:13, 13 May 2009
You are encouraged to solve this task according to the task description, using any language you may know.
Find an approximating polynom of known degree for a given data.
Example: For input data:
x = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10}; y = {1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321};
The approximating polynom is:
3 x2 + 2 x + 1
Here, the polynom's coefficients are (3, 2, 1).
This task is intended as a subtask for Measure relative performance of sorting algorithms implementations.
Ada
<lang ada> with Ada.Numerics.Real_Arrays; use Ada.Numerics.Real_Arrays;
function Fit (X, Y : Real_Vector; N : Positive) return Real_Vector is
A : Real_Matrix (0..N, X'Range); -- The plane
begin
for I in A'Range (2) loop for J in A'Range (1) loop A (J, I) := X (I)**J; end loop; end loop; return Solve (A * Transpose (A), A * Y);
end Fit; </lang> The function Fit implements least squares approximation of a function defined in the points as specified by the arrays xi and yi. The basis φj is xj, j=0,1,..,N. The implementation is straightforward. First the plane matrix A is created. Aji=φj(xi). Then the linear problem AATc=Ay is solved. The result cj are the coefficients. Constraint_Error is propagated when dimensions of X and Y differ or else when the problem is ill-defined.
Example
<lang ada> with Fit; with Ada.Float_Text_IO; use Ada.Float_Text_IO;
procedure Fitting is
C : constant Real_Vector := Fit ( (0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0), (1.0, 6.0, 17.0, 34.0, 57.0, 86.0, 121.0, 162.0, 209.0, 262.0, 321.0), 2 );
begin
Put (C (0), Aft => 3, Exp => 0); Put (C (1), Aft => 3, Exp => 0); Put (C (2), Aft => 3, Exp => 0);
end Fitting; </lang> Sample output:
1.000 2.000 3.000
ALGOL 68
<lang algol>MODE FIELD = REAL;
MODE
VEC = [0]FIELD, MAT = [0,0]FIELD;
PROC VOID raise index error := VOID: (
print(("stop", new line)); stop
);
COMMENT from http://rosettacode.org/wiki/Matrix_Transpose#ALGOL_68 END COMMENT OP ZIP = ([,]FIELD in)[,]FIELD:(
[2 LWB in:2 UPB in,1 LWB in:1UPB in]FIELD out; FOR i FROM LWB in TO UPB in DO out[,i]:=in[i,] OD; out
);
COMMENT http://rosettacode.org/wiki/Matrix_multiplication#ALGOL_68 END COMMENT OP * = (VEC a,b)FIELD: ( # basically the dot product #
FIELD result:=0; IF LWB a/=LWB b OR UPB a/=UPB b THEN raise index error FI; FOR i FROM LWB a TO UPB a DO result+:= a[i]*b[i] OD; result );
OP * = (VEC a, MAT b)VEC: ( # overload vector times matrix #
[2 LWB b:2 UPB b]FIELD result; IF LWB a/=LWB b OR UPB a/=UPB b THEN raise index error FI; FOR j FROM 2 LWB b TO 2 UPB b DO result[j]:=a*b[,j] OD; result );
OP * = (MAT a, b)MAT: ( # overload matrix times matrix #
[LWB a:UPB a, 2 LWB b:2 UPB b]FIELD result; IF 2 LWB a/=LWB b OR 2 UPB a/=UPB b THEN raise index error FI; FOR k FROM LWB result TO UPB result DO result[k,]:=a[k,]*b OD; result );
COMMENT from http://rosettacode.org/wiki/Pyramid_of_numbers#ALGOL_68 END COMMENT OP / = (VEC a, MAT b)VEC: ( # vector division #
[LWB a:UPB a,1]FIELD transpose a; transpose a[,1]:=a; (transpose a/b)[,1]
);
OP / = (MAT a, MAT b)MAT:( # matrix division #
[LWB b:UPB b]INT p ; INT sign; [,]FIELD lu = lu decomp(b, p, sign); [LWB a:UPB a, 2 LWB a:2 UPB a]FIELD out; FOR col FROM 2 LWB a TO 2 UPB a DO out[,col] := lu solve(b, lu, p, a[,col]) [@LWB out[,col]] OD; out
);
FORMAT int repr = $g(0)$,
real repr = $g(-7,4)$;
PROC fit = (VEC x, y, INT order)VEC: BEGIN
[0:order, LWB x:UPB x]FIELD a; # the plane # FOR i FROM 2 LWB a TO 2 UPB a DO FOR j FROM LWB a TO UPB a DO a [j, i] := x [i]**j OD OD; ( y * ZIP a ) / ( a * ZIP a )
END # fit #;
PROC print polynomial = (VEC x)VOID: (
BOOL empty := TRUE; FOR i FROM UPB x BY -1 TO LWB x DO IF x[i] NE 0 THEN IF x[i] > 0 AND NOT empty THEN print ("+") FI; empty := FALSE; IF x[i] NE 1 OR i=0 THEN IF ENTIER x[i] = x[i] THEN printf((int repr, x[i])) ELSE printf((real repr, x[i])) FI FI; CASE i+1 IN SKIP,print(("x")) OUT printf(($"x**"g(0)$,i)) ESAC FI OD; IF empty THEN print("0") FI; print(new line)
);
fitting: BEGIN
VEC c = fit ( (0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0), (1.0, 6.0, 17.0, 34.0, 57.0, 86.0, 121.0, 162.0, 209.0, 262.0, 321.0), 2 ); print polynomial(c); VEC d = fit ( (0, 1, 2, 3, 4, 5, 6, 7, 8, 9), (2.7, 2.8, 31.4, 38.1, 58.0, 76.2, 100.5, 130.0, 149.3, 180.0), 2 ); print polynomial(d)
END # fitting #</lang> output:
3x**2+2x+1 1.0848x**2+10.3552x-0.6164
C
Include file (to make the code reusable easily) named polifitgsl.h <lang c>#ifndef _POLIFITGSL_H
- define _POLIFITGSL_H
- include <gsl/gsl_multifit.h>
- include <stdbool.h>
- include <math.h>
bool polynomialfit(int obs, int degree, double *dx, double *dy, double *store); /* n, p */
- endif</lang>
Implementation (the examples here helped alot to code this quickly): <lang c>#include "polifitgsl.h"
bool polynomialfit(int obs, int degree, double *dx, double *dy, double *store) /* n, p */ {
gsl_multifit_linear_workspace *ws; gsl_matrix *cov, *X; gsl_vector *y, *c; double chisq;
int i, j;
X = gsl_matrix_alloc(obs, degree); y = gsl_vector_alloc(obs); c = gsl_vector_alloc(degree); cov = gsl_matrix_alloc(degree, degree);
for(i=0; i < obs; i++) { gsl_matrix_set(X, i, 0, 1.0); for(j=1; j < degree; j++) { gsl_matrix_set(X, i, j, pow(dx[i], j)); } gsl_vector_set(y, i, dy[i]); }
ws = gsl_multifit_linear_alloc(obs, degree); gsl_multifit_linear(X, y, c, cov, &chisq, ws);
/* store result ... */ for(i=0; i < degree; i++) { store[i] = gsl_vector_get(c, i); }
gsl_multifit_linear_free(ws); gsl_matrix_free(X); gsl_matrix_free(cov); gsl_vector_free(y); gsl_vector_free(c); return true; /* we do not "analyse" the result (cov matrix mainly)
to know if the fit is "good" */ }</lang> Testing: <lang c>#include <stdio.h>
- include "polifitgsl.h"
- define NP 11
double x[] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10}; double y[] = {1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321};
- define DEGREE 3
double coeff[DEGREE];
int main() {
int i;
polynomialfit(NP, DEGREE, x, y, coeff); for(i=0; i < DEGREE; i++) { printf("%lf\n", coeff[i]); } return 0;
}</lang> Output of the test:
1.000000 2.000000 3.000000
Fortran
<lang fortran>
module fitting contains
function polyfit(vx, vy, d) implicit none integer, intent(in) :: d integer, parameter :: dp = selected_real_kind(15, 307) real(dp), dimension(d+1) :: polyfit real(dp), dimension(:), intent(in) :: vx, vy real(dp), dimension(:,:), allocatable :: X real(dp), dimension(:,:), allocatable :: XT real(dp), dimension(:,:), allocatable :: XTX integer :: i, j integer :: n, lda, lwork integer :: info integer, dimension(:), allocatable :: ipiv real(dp), dimension(:), allocatable :: work n = d+1 lda = n lwork = n allocate(ipiv(n)) allocate(work(lwork)) allocate(XT(n, size(vx))) allocate(X(size(vx), n)) allocate(XTX(n, n)) ! prepare the matrix do i = 0, d do j = 1, size(vx) X(j, i+1) = vx(j)**i end do end do XT = transpose(X) XTX = matmul(XT, X) ! calls to LAPACK subs DGETRF and DGETRI call DGETRF(n, n, XTX, lda, ipiv, info) if ( info /= 0 ) then print *, "problem" return end if call DGETRI(n, XTX, lda, ipiv, work, lwork, info) if ( info /= 0 ) then print *, "problem" return end if polyfit = matmul( matmul(XTX, XT), vy) deallocate(ipiv) deallocate(work) deallocate(X) deallocate(XT) deallocate(XTX) end function end module
</lang>
Example
<lang fortran>
program PolynomalFitting use fitting implicit none ! let us test it integer, parameter :: degree = 2 integer, parameter :: dp = selected_real_kind(15, 307) integer :: i real(dp), dimension(11) :: x = (/ (i,i=0,10) /) real(dp), dimension(11) :: y = (/ 1, 6, 17, 34, & 57, 86, 121, 162, & 209, 262, 321 /) real(dp), dimension(degree+1) :: a a = polyfit(x, y, degree) write (*, '(F9.4)') a
end program
</lang>
Output (lower powers first, so this seems the opposite of the Python output):
1.0000 2.0000 3.0000
Octave
<lang octave>x = [0:10]; y = [1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321]; coeffs = polyfit(x, y, 2)</lang>
Python
<lang python> >>> x = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10] >>> y = [1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321] >>> coeffs = numpy.polyfit(x,y,deg=2) >>> coeffs array([ 3., 2., 1.]) </lang> Substitute back received coefficients. <lang python> >>> yf = numpy.polyval(numpy.poly1d(coeffs), x) >>> yf array([ 1., 6., 17., 34., 57., 86., 121., 162., 209., 262., 321.]) </lang> Find max absolute error. <lang python> >>> '%.1g' % max(y-yf) '1e-013' </lang>
Example
For input arrays `x' and `y': <lang python> >>> x = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9] >>> y = [2.7, 2.8, 31.4, 38.1, 58.0, 76.2, 100.5, 130.0, 149.3, 180.0] </lang>
<lang python> >>> p = numpy.poly1d(numpy.polyfit(x, y, deg=2), variable='N') >>> print p
2
1.085 N + 10.36 N - 0.6164 </lang> Thus we confirm once more that for already sorted sequences the considered quick sort implementation has quadratic dependence on sequence length (see Example section for Python language on Query Performance page).