Algoritmi

L'algoritmo dei minimi quadrati

Approssimazione di una serie di dati mediante una serie di funzioni con errore quadratico minimo.

Il metodo dei minimi quadrati consente di approssimare mediante una serie di funzioni una serie di dati con errore quadratico minimo; quando la serie di funzioni considerata è costituita da polinomi di grado non superiore al primo si afferma che la curva di regressione è lineare.
Per comprendere meglio l'algoritmo che vogliamo implementare analizziamo il grafico di un risultato:

Nella figura precedente si notano i punti rappresentativi dei dati e la retta dell'approssimazione ottenuta con il metodo dei minimi quadrati.

Consideriamo il metodo di costruzione della curva di regressione avente come base dello spazio di funzioni l’insieme {1, x}, tuttavia il metodo dei minimi quadrati consente di costruire una curva di regressione per combinazione lineare di uno spazio di funzioni qualsiasi.
Sappiamo che per due punti passa una ed una sola retta, considerando una serie di n punti con n>2, non è possibile costruire una retta di interpolazione a meno che volutamente o fortuitamente la serie di n punti abbia andamento lineare, quello che potremo però fare è costruire la retta che meglio approssima la serie di n punti.

In termini più formali, avendo a disposizione la serie di n punti:

vogliamo determinare la retta y=mx+q tale che la somma degli scarti quadratici dai punti della serie sia minima.

Lo scarto quadratico è definito pari a (y_i-y(x_i))², ovvero la differenza tra il valore reale della serie e quello stimato dalla retta di regressione.

Il residuo è definito come la somma degli scarti quadratici medi:

    [1]

ed y(x_i) = mx_i+q  è  il valore stimato dalla retta di regressione lineare nel punto x_i.

Dall’analisi matematica, si ricava che la retta per cui il residuo è minimo è fornita dal minimo locale della equazione di residuo [1] rispetto alle due variabili m e q.

La tecnica di ricerca di un minimo locale implica la risoluzione del sistema di equazioni ottenuto uguagliando a zero le derivate parziali della [1]:

                    [2]

Si ottiene così un sistema lineare di due equazioni e due incognite, nelle variabili m e q:

       [3]

si osservi che le uniche variabili sono m e q, mentre tutti gli altri termini sono numeri calcolabili dai dati in ingresso, ovvero:

n: numero dei dati in ingresso

: somma delle ascisse dei punti noti
: somma dei quadrati delle ascisse dei punti noti

:somma delle ordinate dei punti noti
:somma dei prodotti delle ascisse per le ordinate dei punti noti

La soluzione del sistema di equazioni [3] è, quindi:

             [4]

Sostituendo alla [4] i valori del caso in esame siamo in grado di scrivere la retta di regressione

y = mx+q

Analizziamo infine un risultato esplicativo:

Di seguito è riportata una routine in fortran 90 per l'applicazione del metodo dei minimi quadrati

!*****************************************************************
!*         LEAST SQUARES POLYNOMIAL FITTING PROCEDURE            *
!* ------------------------------------------------------------- *
!* This program least squares fits a polynomial to input data.   *
!* forsythe orthogonal polynomials are used in the fitting.      *
!* The number of data points is n.                               *
!* The data is input to the subroutine in x[i], y[i] pairs.      *
!* The coefficients are returned in c[i],                        *
!* the smoothed data is returned in v[i],                        *
!* the order of the fit is specified by m.                       *
!* The standard deviation of the fit is returned in d.           *
!* There are two options available by use of the parameter e:    *
!*  1. if e = 0, the fit is to order m,                          *
!*  2. if e > 0, the order of fit increases towards m, but will  *
!*     stop if the relative standard deviation does not decrease *
!*     by more than e between successive fits.                   *
!* The order of the fit then obtained is l.                      *
!*****************************************************************
Subroutine LS_POLY(m,e1,n,l,x,y,c,dd)
  parameter(SIZE=25)
  !Labels: 10,15,20,30,50
  real*8 x(SIZE),y(SIZE),v(SIZE),a(SIZE),b(SIZE)
  real*8 c(0:SIZE),d(SIZE),c2(SIZE),e(SIZE),f(SIZE)
  integer i,l,l2,m,n,n1
  real*8 a1,a2,b1,b2,c1,dd,d1,e1,f1,f2,v1,v2,w
  n1 = m + 1; l=0
  v1 = 1.d7
  ! Initialize the arrays
  do i=1, n1
    a(i) = 0.d0; b(i) = 0.d0; f(i) = 0.d0
  end do
  do i=1, n
    v(i) = 0.d0; d(i) = 0.d0
  end do
  d1 = dsqrt(dfloat(n)); w = d1;
  do i=1, n
    e(i) = 1.d0 / w
  end do
  f1 = d1; a1 = 0.d0
  do i=1, n
    a1 = a1 + x(i) * e(i) * e(i)
  end do
  c1 = 0.d0
  do i=1, n
    c1 = c1 + y(i) * e(i)
  end do
  b(1) = 1.d0 / f1; f(1) = b(1) * c1
  do i=1, n
    v(i) = v(i) + e(i) * c1
  end do
  m = 1
! Save latest results
10 do i=1, l
    c2(i) = c(i)
  end do
  l2 = l; v2 = v1; f2 = f1; a2 = a1; f1 = 0.d0
  do i=1, n
    b1 = e(i)
    e(i) = (x(i) - a2) * e(i) - f2 * d(i)
    d(i) = b1
    f1 = f1 + e(i) * e(i)
  end do
  f1 = dsqrt(f1)
  do i=1, n
    e(i) = e(i) / f1
  end do
  a1 = 0.d0
  do i=1, n  
    a1 = a1 + x(i) * e(i) * e(i)
  end do
  c1 = 0.d0
  do i=1, n  
    c1 = c1 + e(i) * y(i)
  end do	 
  m = m + 1; i = 0
15 l = m - i; b2 = b(l); d1 = 0.d0
  if (l > 1)  d1 = b(l - 1)
  d1 = d1 - a2 * b(l) - f2 * a(l)
  b(l) = d1 / f1; a(l) = b2; i = i + 1
  if (i.ne.m) goto 15
  do i=1, n
    v(i) = v(i) + e(i) * c1 
  end do
  do i=1, n1
    f(i) = f(i) + b(i) * c1
    c(i) = f(i)
  end do
  vv = 0.d0
  do i=1, n
    vv = vv + (v(i) - y(i)) * (v(i) - y(i))
  end do	 
  !Note the division is by the number of degrees of freedom
  vv = dsqrt(vv / dfloat(n - l - 1)); l = m
  if (e1.eq.0.d0) goto 20
  !Test for minimal improvement
  if (dabs(v1 - vv) / vv < e1) goto 50
  !if error is larger, quit
  if (e1 * vv > e1 * v1) goto 50
  v1 = vv
20 if (m.eq.n1) goto 30
  goto 10
!Shift the c[i] down, so c(0) is the constant term
30 do i=1, l  
    c(i - 1) = c(i)
  end do
  c(l) = 0.d0
  ! l is the order of the polynomial fitted
  l = l - 1; dd = vv
  return
! Aborted sequence, recover last values
50 l = l2; vv = v2
  do i=1, l  
    c(i) = c2(i)
  end do
  goto 30
end

• » Invia tramite EMAIL
• » Versione per la STAMPA
• » Le vostre opinioni

I link correlati all'argomento

• L'algoritmo di Eulero
• L'algoritmo di eliminazione di Gauss
• L'algoritmo della Secante
• L'algoritmo della bisezione