Esercizio 2

Utilizzando i minimi quadrati stimare i valori $A$ e $B$ della retta migliore $y=A+Bx$ rappresentata dalle nostra doppia lista di $N$ dati (le coppie ($x_i,y_i$)). Ricordiamo che le equazioni che dobbiamo calcolare sono

$$\begin{eqnarray*} D=N(\sum_{i=1}^{N}x_i^2) -(\sum_{i=1}^{N}x_i)^2\\ A=\frac{(\s... ...sum_{i=1}^{N}x_i y_i)-(\sum_{i=1}^{N}x_i)(\sum_{i=1}^{N}y_i)}{D} \end{eqnarray*}$$

Risoluzione: innanzi tutto, abbiamo le equazioni precedenti come forma matematica della soluzione del nostro algoritmo. Il primo punto da affrontare sarà quindi il tipo di struttura di dati che ci serve. Dalle equazioni precedenti, vediamo che abbiamo bisogno di due insiemi di $N$ dati, le nostre $x$ e le corrispondenti $y$. Una possibile soluzione è quindi quella di utilizzare due liste, x=[], y=[] che ci servono per utilizzare all'interno del nostro programma i valori delle $x$ e $y$.

Seguendo un approccio incrementale allora prima scriviamo un programma che legga i dati, li immagazzini nelle nostre strutture dati (liste) e poi stampi il contenuto di ciò che ha letto. In questo modo, ci curiamo che la parte di input/output sia corretta e non la toccheremo più. In secondo luogo, possiamo notare che le equazioni precedenti si possono scomporre in

• somma degli elementi della lista ($\sum x_i$ e $\sum y_i$)
• somma dei prodotti dei corrispondenti elementi di due liste (o della medesima) ( $\sum x_i^2=\sum x_ix_i$ e $\sum x_iy_i$)

che possiamo identificare come sumX, sumY e sumXX e sumXY. Da cui il nostro primo programma che legge soltanto la lista e ne restituisce il valore può essere

if __name__ == '__main__':
   numdata=input('number of data = ')
   i = 0
   x=[]
   y=[]
   while i < numdata :
       # input legge anche una lista di valori separati da ,
       d1,d2=input('give me the x,y value = ') 
       x.append(d1)
       y.append(d2)
       i =i+1
   N=numdata
   A = 0.0
   B = 0.0 
   '''
   # equations
   D=N*sumXX - (sumX)**2.0
   A=(sumXX*sumY-sumX*sumXY)/D
   B=(N*sumXY-sumX*sumY)/D
   '''
   print "data",x,y
   print "B x + A => B =",B, "A =", A

A questo punto possiamo fare una importante considerazione. Siccome abbiamo bisogno della somma degli elemeti della lista $x$ e $y$, allora è sufficiente scrivere una sola funzione da utilizzare per entrambe le liste. Possiamo quindi in modo incrementale, scrivere la nostra funzione somma e verificare che faccia quello per cui l'abbiamo definita. Da cui scriviamo la funzione

def sum(v):
    ''' sum(v) -> sum over the element v[i] '''
    s=0
    for e in v:
        s = s + e
    return float(s)

e nel programma principale aggiungiamo le istruzioni

    sumX=sum(x)
    sumY=sum(y)
    print 'sumX=',sumX,'sumY='sumY # to test the function

Analogamente al caso precedente, possiamo notare che anche per sumXX e sumXY è necessaria una sola funzione, per esempio

def sum2(v1,v2):
    ''' sum2(v1,v2) -> sum over the v1[i]*v2[i] 
    '''
    d=len(v1)
    s=0
    for i in range(d):
        s = s + v1[i]*v2[i]
    return float(s)

di cui possiamo verificare il comportamento, inserendo nel programma le seguenti istruzioni

   sumXX=sum2(x,x)
   sumXY=sum2(x,y)
   print 'sumXX=',sumXX,'sumXY='sumXY # to test the function

Una volta verificato che i singoli pezzi del nostro programma singolarmente fanno il loro dovere, allora possiamo vedere se le singole parti funzionano bene assieme, per cui il nostro programma finale sarà

def sum(v):
    ''' sum(v) -> sum over the element v[i] '''
    s=0
    for e in v:
        s = s + e
    return float(s)

def sum2(v1,v2):
    ''' sum2(v1,v2) -> sum over the v1[i]*v2[i] 
    '''
    d=len(v1)
    s=0
    for i in range(d):
        s = s + v1[i]*v2[i]
    return float(s)

if __name__ == '__main__':
   numdata=input('number of data = ')
   i = 0
   x=[]
   y=[]
   while i < numdata :
       d1,d2=input('give me the x,y value = ')
       x.append(d1)
       y.append(d2)
       i =i+1
   N=numdata
   sumX=sum(x)
   sumY=sum(y)
   sumXX=sum2(x,x)
   sumXY=sum2(x,y)
   # equations
   D=N*sumXX - (sumX)**2.0
   A=(sumXX*sumY-sumX*sumXY)/D
   B=(N*sumXY-sumX*sumY)/D
   print "data",x,y
   print "B x + A => B =",B, "A =", A