Massimizzazione di funzioni

Nella sezione precedente abbiamo visto che per ottenere i valori che massimizzano la likelihood deriviamo rispetto ai paramertri della probabilità che vogliamo massimizzare, e eguagliamo a zero il il risultato. A volte è però difficile (se non impossibile) ottenere una semplice soluzione analitica come abbiamo visto per i casi precedenti, nonostante ciò siamo in grado di trovare soluzioni numeriche al nostro problema.

Esistono diversi algoritmi iterativi per otenere il risultato, il più semplice (e quello che vedremo) è basato sulla discesa/ascesa lungo il gradiente. Se infatti ricordiamo che il gradiente misura la variazione della pendenza rispetto alle varie coordinate siamo in grado di aggiornare i nostri parametri secondo la semplice equazione

$$\begin{eqnarray*} p_n = p_{n-1} + s\frac{\partial f(p,x)}{\partial p} \end{eqnarray*}$$

dove $p$ è il nostro parametro di interesse e $s$ è una variabile arbitraria che ci dice quanto grande deve essere il nostro spostamento per il parametro $p$ allo step $n$ nella direzione del gradiente.

Per vedere come questo funziona in pratica, riprendiamo l'esempio della gaussiana di equazione . Abbiamo già visto come la soluzione analitica della massimizzazione della likelihood nel caso di un insieme di $n$ dati fornisca proprio la definizione del valor medio e deviazione standard alla quale siamo abituati. Possiamo quindi costruire la nostra funzione logaritmo della gaussiana (lf(x,parameters)) e le sue derivate parziali (partiali_mu, partial_sigma). Per cui segue che una possibile implementazione della equazione precedente diviene la deltamax, ossia

def lf(x,parameters):
    ''' ln(gaussian) in x with parameter=mu, sigma   
    '''
    mu,sigma=parameters
    return -math.log(sigma) -((x-mu)**2)/(2*sigma*sigma) -0.5*math.log(2*math.pi)

def partial_mu(x,parameters):
    ''' partial derivative of the gaussian with respect to mu in x '''
    mu,sigma=parameters
    return (x-mu)/(sigma*sigma)

def partial_sigma(x,parameters):
    ''' partial derivative of the gaussian with respect to sigma in x'''
    mu,sigma=parameters
    return ( ((x-mu)*(x-mu))/(sigma*sigma) -1 )/sigma

def deltamax(f,partials,parameters,vars,delta=0.001,tol=0.0001,maxcyc=300):
    ''' maximize using explicit partial derivativesi
        f => function to maximize
        partials => list of parial derivative funcions
        parameters => function parameters
        vars => list of x values
        delta => modification step
        tol => tolerance to which stop
        maxcyc => maximum number of allowed cycles 
    '''
    assert len(partials) == len(parameters)
    oldpar=parameters[:]
    newpar=parameters[:]
    valpartials=[0.0]*len(partials)
    err=tol+100
    c=0
    while tol < err and c < maxcyc:
        err=0.0
        for i in range(len(parameters)):
            valpartials[i]=0.0
            for x in vars:
                 valpartials[i] += partials[i](x,oldpar)
            newpar[i]=oldpar[i]+delta*valpartials[i]
            err=err+abs(newpar[i]-oldpar[i])
        oldpar=newpar[:]
        c=c+1
    return oldpar

Se invecie non siamo in grado di calcolare le derivate, oppure non ne abbiamo voglia per nulla, allora possiamo come in precedenza fare uso della definizione di derivata per calcolare le derivate della nostra funzione senza dover perdere tempo con le tavole delle derivate. In questo caso il codice diviene ancora più semplice e più generale (anche se bisogna stare molto attenti allo step di integrazione che può essere necessario sia diveso per ogni parametro) come possiamo osservare

def numderivative(f,x,parameters,idx_par,step=0.0001):
    ''' 
       numeric derivative of f with respect to parameter 
       whose index is idx_par 
    '''
    nextpar=parameters[:]
    nextpar[idx_par]+=step
    return (f(x,nextpar) - f(x,parameters))/step


def deltaNumMax(f,parameters,vars,delta=0.001,tol=0.0001,maxcyc=300):
    ''' maximize using numeric derivatives
        f => function to maximize
        parameters => function parameters
        vars => list of x values
        delta => modification step
        tol => tolerance to which stop
        maxcyc => maximum number of allowed cycles 
    '''
    oldpar=parameters[:]
    newpar=parameters[:]
    valpartials=[0.0]*len(parameters)
    err=tol+100
    c=0
    while tol < err and c < maxcyc:
        err=0.0
        for i in range(len(parameters)):
            valpartials[i]=0.0
            for x in vars:
                 valpartials[i] += numderivative(f,x,oldpar,i) 
            newpar[i]=oldpar[i]+delta*valpartials[i]
            err=err+abs(newpar[i]-oldpar[i])
        oldpar=newpar[:]
        c=c+1
    return oldpar

Infine essendo noi scettici, vogliamo assolutamente verificare che questo codice corrisponda a quanto abbiamo affermato, per cui possiamo vedere se è vero che otteniamo lo stesso risultaro del valor medio e deviazione standard. Per farlo potremo costruire il semplice codice

if __name__=='__main__':
     import random
     v=[]
     mu=0.1
     sigma=2.0
     m=1000
     param=[mu, sigma]
     for i in range(m):
        v.append(random.normalvariate(0.1,2.0))
     print "true equations",
     print average(v),sig2(v),math.sqrt(sig2(v))
     print "Explicit derivatives", 
     print deltamaxi(lf,[partial_mu,partial_sigma],[0.0,10.0],v)
     print "Numeric derivatives",
     print deltaNumMax(lf,[0.0,10.0],v)

e valutare se forniscono il medesimo risultato.

La cosa importante è che ora abbiamo uno strumento indipendente dalla complessità della funzione utilizzata che ci permette di trovare soluzioni che massimizzano la nostra probabilià.

Molto importante: dobbiamo però ricordare che nel caso generale quello che noi otteniamo è un massimo locale, in quanto partiamo da dei precisi valori dei nostri parametri ([0.0,10.0] nel nostro caso) e che quindi questo procedimento non garantisce che la soluzione sia globale.