Polynomial regression: Difference between revisions

From Rosetta Code
Content added Content deleted
Line 60: Line 60:
<lang algol>MODE FIELD = REAL;
<lang algol>MODE FIELD = REAL;


MODE
MODE
VEC = [0]FIELD,
VEC = [0]FIELD,
MAT = [0,0]FIELD;
MAT = [0,0]FIELD;
Line 73: Line 73:
[2 LWB in:2 UPB in,1 LWB in:1UPB in]FIELD out;
[2 LWB in:2 UPB in,1 LWB in:1UPB in]FIELD out;
FOR i FROM LWB in TO UPB in DO
FOR i FROM LWB in TO UPB in DO
out[,i]:=in[i,]
out[,i]:=in[i,]
OD;
OD;
out
out
);
);


COMMENT http://rosettacode.org/wiki/Matrix_multiplication#ALGOL_68 END COMMENT
COMMENT from http://rosettacode.org/wiki/Matrix_multiplication#ALGOL_68 END COMMENT
OP * = (VEC a,b)FIELD: ( # basically the dot product #
OP * = (VEC a,b)FIELD: ( # basically the dot product #
FIELD result:=0;
FIELD result:=0;
Line 112: Line 112:
[,]FIELD lu = lu decomp(b, p, sign);
[,]FIELD lu = lu decomp(b, p, sign);
[LWB a:UPB a, 2 LWB a:2 UPB a]FIELD out;
[LWB a:UPB a, 2 LWB a:2 UPB a]FIELD out;
FOR col FROM 2 LWB a TO 2 UPB a DO
FOR col FROM 2 LWB a TO 2 UPB a DO
out[,col] := lu solve(b, lu, p, a[,col]) [@LWB out[,col]]
out[,col] := lu solve(b, lu, p, a[,col]) [@LWB out[,col]]
OD;
OD;
Line 172: Line 172:
print polynomial(d)
print polynomial(d)
END # fitting #</lang>
END # fitting #</lang>
Output:
output:
<pre>
<pre>
3x**2+2x+1
3x**2+2x+1
1.0848x**2+10.3552x-0.6164
1.0848x**2+10.3552x-0.6164
</pre>
</pre>

=={{header|C}}==
=={{header|C}}==
{{libheader|libgsl}}
{{libheader|libgsl}}

Revision as of 13:18, 13 May 2009

Task
Polynomial regression
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

Translation of: Ada
Works with: ALGOL 68 version Standard - lu decomp and lu solve are from the GSL library
Works with: ALGOL 68G version Any - tested with release mk15-0.8b.fc9.i386

<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 from 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

Library: libgsl

Include file (to make the code reusable easily) named polifitgsl.h <lang c>#ifndef _POLIFITGSL_H

  1. define _POLIFITGSL_H
  2. include <gsl/gsl_multifit.h>
  3. include <stdbool.h>
  4. include <math.h>

bool polynomialfit(int obs, int degree, double *dx, double *dy, double *store); /* n, p */

  1. 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>

  1. include "polifitgsl.h"
  1. 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};

  1. 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

Library: LAPACK

<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

Library: numpy

<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).

Retrieved from "https://rosettacode.org/wiki/Polynomial_regression?oldid=53486"