come minimizzare una funzione con valori di variabile discreti in scipy

Question

Sto tentando di ottimizzare una funzione di destinazione che ha più variabili di input (tra 24 e 30). Queste variabili sono esempi di tre diverse variabili statistiche e i valori delle funzioni di destinazione sono valori di probabilità del test t. Una funzione di errore rappresenta l'errore (somma dei quadrati delle differenze) tra le probabilità di t-test desiderate e quelle effettive. Posso accettare solo soluzioni in cui l'errore è inferiore a 1e-8, per tutti e tre i t-test.

Stavo usando scipy.optimize.fmin e ha funzionato benissimo. Esistono molte soluzioni in cui la funzione obiettivo è diventata zero.

Il problema è che ho bisogno di trovare una soluzione in cui le variabili sono tra 0 e 10.0, e sono numeri interi o non hanno più di una cifra della parte frazionaria. Esempi di valori validi sono 0 10 3 5.5 6.8. Esempi di valori non validi: -3 2.23 30 o 0.16666667.

Mi capita di sapere che esiste almeno una soluzione, perché i valori obiettivo provengono da dati misurati effettivi. I dati originali sono andati persi e il mio compito è trovarli. Ma non so come. L'utilizzo di tentativi/errori non è un'opzione, poiché ci sono circa 100 valori possibili per ciascuna variabile e, dato il numero di variabili, il numero di casi possibili sarebbe 100 ** 30, che è troppo. L'utilizzo di fmin è ottimo, tuttavia non funziona con valori discreti.

C'è un modo per risolvere questo? Non è un problema se ho bisogno di eseguire un programma per molte ore per trovare una soluzione. Ma ho bisogno di trovare soluzioni per circa 10 valori target entro pochi giorni, e sono fuori da nuove idee.

Ecco un esempio MWE:

import math 
import numpy 
import scipy.optimize 
import scipy.stats 
import sys 

def log(s): 
    sys.stdout.write(str(s)) 
    sys.stdout.flush() 

# List of target T values: TAB, TCA, TCB 
TARGETS = numpy.array([ 
    [0.05456834, 0.01510358, 0.15223353 ], # task 1 to solve 
    [0.15891875, 0.0083665,  0.00040262 ], # task 2 to solve 
]) 
MAX_ERR = 1e-10 # Maximum error in T values 
NMIN,NMAX = 8,10 # Number of samples for T probes. Inclusive. 

def fsq(x, t, n): 
    """Returns the differences between the target and the actual values.""" 
    a,b,c = x[0:n],x[n:2*n],x[2*n:3*n] 
    results = numpy.array([ 
     scipy.stats.ttest_rel(a,b)[1], # ab 
     scipy.stats.ttest_rel(c,a)[1], # ca 
     scipy.stats.ttest_rel(c,b)[1] # cb 
    ]) 
    # Sum of squares of diffs 
    return (results - t) 

def f(x, t, n): 
    """This is the target function that needs to be minimized.""" 
    return (fsq(x,t,n)**2).sum() 

def main(): 
    for tidx,t in enumerate(TARGETS): 
     print "=============================================" 
     print "Target %d/%d"%(tidx+1,len(TARGETS)) 
     for n in range(NMIN,NMAX+1): 
      log(" => n=%s "%n) 
      successful = False 
      tries = 0 
      factor = 0.1 
      while not successful: 
       x0 = numpy.random.random(3*n) * factor 
       x = scipy.optimize.fmin(f,x0, [t,n], xtol=MAX_ERR, ftol=MAX_ERR) 
       diffs = fsq(x,t,n) 
       successful = (numpy.abs(diffs)<MAX_ERR).all() 
       if successful: 
        log(" OK, error=[%s,%s,%s]\n"%(diffs[0],diffs[1],diffs[2])) 
        print " SOLUTION FOUND " 
        print x 
       else: 
        tries += 1 
        log(" FAILED, tries=%d\n"%tries) 
        print diffs 
        factor += 0.1 
        if tries>5: 
         print "!!!!!!!!!!!! GIVING UP !!!!!!!!!!!" 
         break 
if __name__ == "__main__": 
    main() 

Answer 1

Cosa si sta tentando di fare (se ho capito la configurazione) si chiama programmazione intera ed è NP-hard; http://en.wikipedia.org/wiki/Integer_programming. Mi rendo conto che non stai cercando soluzioni intere, ma se moltiplichi tutti i tuoi input per 10 e dividi la funzione target per 100, ottieni un problema equivalente dove gli input sono tutti interi. Il punto è che i tuoi input sono discreti.

La funzione obiettivo su cui si sta lavorando è una funzione convessa, quadratica, e ci sono buoni algoritmi di ottimizzazione vincolati che la risolveranno rapidamente per gli input a valori reali nell'intervallo [0, 10]. Da questo puoi provare a arrotondare o controllare tutti i punti accettabili nelle vicinanze, ma ce ne sono 2^n, dove n è il numero di input. Anche se lo fai, la soluzione ottimale non è garantita per essere uno di questi punti.

Ci sono algoritmi di approssimazione per problemi di programmazione di interi e potreste scoprire che a volte uno di essi funziona abbastanza bene per portarvi a un punto ottimale. Ci sono un elenco di cose che potresti provare nell'articolo di Wikipedia che ho citato, ma non so che sarai felice nel cercare di risolvere questo problema.