TI/Programowanie dla Fizyków Medycznych/Optymalizacja
Spis treści
Optymalizacja jednowymiarowa
Omawianie zagadnienia optymalizacji rozpocznijmy od prostego przykładu. Zdefiniujmy pewną funkcję i zobaczmy jak wygląda jej wykres.
import numpy as np
import pylab as py
licznikTestowej=0
def testowa(x):
global licznikTestowej
licznikTestowej+=1
return 1/x+np.exp(x)
xtest=np.arange(0.2,2,0.01)
ytest=[testowa(x) for x in xtest]
py.plot(xtest,ytest)
py.show()
Na rozważanym przedziale [0.2,2] powyższa funkcja ma tylko jedno ekstremum lokalne. Taką funkcję nazywamy unimodalną. Zmienna licznikTestowej umożliwi nam zliczanie wywołań funkcji testowej przez analizowane procedury. Zagadnienie, którym teraz będziemy się zajmować to problem numerycznego znajdowania takiego ekstremum. Jak w każdym problemie numerycznym ekstremum szukać będziemy zakładając pewną dokładność otrzymanego wyniku, którą oznaczmy xtol. Na wstępie przyjmijmy, że poszukujemy ekstremum z dokładnością xtol=0.01. Najprostszą metodą będzie policzenie wartości funkcji dla wszystkich wartości x z podanego przedziału co xtol. Jest to metoda siłowa i wielokrotnie licząca wartość funkcji. Jej kod możemy znaleźć poniżej.
def bruteForce(func,xmin,xmax,args=(),xtol=0.01):
xlist=np.arange(xmin,xmax,xtol)
ylist=[func(x,*args) for x in xlist]
return xlist[ylist.index(max(ylist))]
Innym, znacznie efektywniejszym sposobem znajdowania minimum może być następująca procedura rekurencyjna:
- podzielmy przedział [xmin,xmax] na 3 równe część: [xmin,xL],[xL,xR] oraz [xR,xmax]
- jeżeli wartość funkcji w xL jest mniejsza od wartości funkcji w xR to powtórz procedurę dla przedziału [xmin,xR]. W przeciwnym przypadku powtórz procedurę dla przedziału [xL,xmax].
- zakończ działanie gdy badany przedział jest krótszy niż xtol
Przykładowa implementacja tej metody wygląda następująco
def twoMidPointsR(func,xmin,xmax,args=(),xtol=0.01):
if xmax-xmin<xtol: return 0.5*(xmax+xmin)
xL=xmin+(xmax-xmin)/3.0
xR=xmax-(xmax-xmin)/3.0
fxL=func(xL,*args)
fxR=func(xR,*args)
if fxL>fxR:
return twoMidPointsR(func,xmin,xR,args,xtol)
else:
return twoMidPointsR(func,xL,xmax,args,xtol)
Tą metodą możemy już szukać minimum z dowolną dokładnością, co nie będzie skutkowało znacznie większym czasem obliczeń. Zauważmy jednak, że w każdej iteracji wartość minimalizowanej funkcji liczona jest dwukrotnie. Jeżeli mamy do czynienia z funkcją, dla której obliczenie pojedynczej wartości jest bardzo czasochłonne to warto się zastanowić czy nie można tego wyniku poprawić. Zauważmy, że dla działania tej metody wcale nie jest konieczne dzielenie badanego odcinka dokładnie na trzy równe części. Można dokonać tego podziału w zupełnie innej proporcji. Warto tak dobrać punkty xL i xR, aby xR pokrywał się z xL (lub xL z xR) w kolejnym kroku iteracji. Jeżeli dodatkowo stworzymy zmienne przechowujące wcześniej liczone wartości funkcji to uda nam się ograniczyć liczbę wywołań funkcji o połowę. Opisana metoda to tak zwana metoda złotego podziału. Przykładowa implementacja wygląda następująco
def GoldenRatioRearch(func,xmin,xmax,args=(),xtol=0.01):
golden=0.5*(np.sqrt(5.0)-1.0)
xL=xmax-golden*(xmax-xmin)
xR=xmin+golden*(xmax-xmin)
fxL=func(xL,*args)
fxR=func(xR,*args)
while xmax-xmin>xtol:
if fxL<fxR:
xmin=xmin
xmax=xR
xR=xL
fxR=fxL
xL=xmax-golden*(xmax-xmin)
fxL=func(xL,*args)
else:
xmin=xL
xmax=xmax
xL=xR
fxL=fxR
xR=xmin+golden*(xmax-xmin)
fxR=func(xR,*args)
return 0.5*(xmax+xmin)
Napisanej metody optymalizacji możemy użyć w celu dopasowania funkcji do danych empirycznych metodą najmniejszych kwadratów. Przykład takiego zastosowania poniżej.
#suma kwadratów
def squares(a,func,xlist,ylist):
return sum([(func(xlist[i],*a)-ylist[i])**2 for i in range(len(xlist))])
#funkcja liniowa
def liniowa(x,a): return x*a
#generujemy przykladowe xlist
xlist=np.arange(0,1,0.001)
#generujemy wartosci funkcji z szumem
ylist=[liniowa(x,1.23)+0.000001*np.random.randn() for x in xlist]
#najlepsze dopasowanie metoda golden ration
print GoldenRatioRearch(squares,0,10,args=(liniowa,xlist,ylist),xtol=0.01)
Optymalizacja wielowymiarowa
Przejście od optymalizacji jedno- do wielowymiarowej fundamentalnie komplikuje problem. Pierwszym problemem jest istnienie tak zwanych punktów siodłowych. Nie istnieją zatem metody które zawsze znajdują szukane minimum, nawet jezeli wiadomo że takie istnieje. Najpopularniejszą metodą jest downhill symplex lub inaczej metoda Neldera-Meada. Z powodu znaczengo stopnia komplikacji nie ędziemy jej samodzielnie implementować, a jedynie poslużymy sie implementacją z biblioteki scipy.optimize. Dziki niej możemy np. dopasowywać funkcję z więcej niż jednym parametrem do danych eksperymentalnych.
import scipy.optimize as so
#suma kwadratów
def squares(a,func,xlist,ylist):
return sum([(func(xlist[i],*a)-ylist[i])**2 for i in range(len(xlist))])
#funkcja liniowa
def liniowa(x,a,b): return x*a+b
#generujemy przykladowe xlist
xlist=np.arange(0,1,0.001)
#generujemy wartosci funkcji z szumem
ylist=[liniowa(x,1.23,-0.73)+0.000001*np.random.randn() for x in xlist]
#najlepsze dopasowanie metoda golden ration
print so.fmin(squares,(1,0),args=(kwadratowa,xlist,ylist))
Oczywiście dopasowywanie możemy przeprowadzać nie tylko metodą najmniejszych kwadratów.
Zadanie - rozklad Cauchy'ego
Wylosuj 1000 liczb z rozkładu Cauchy'ego z parametrami loc=1.23 i scale=2.0. Do wylosowanych danych dopasuj rozkład Cauche'ego trzema metodami
- METODA 1 - stwórz histogram otrzymanych wartości, znormalizuj go i metodą najmniejszych kwadratów dopsauj gęstość rozkladu do histogramu
- METODA 2 - dopasuj gęstość rozkładu do wylosowanych danych metodą największej wiarygodności
- METODA 3 - z wylosowanych danych stwórz dystrybuante empiryczną. Metodą najmniejszych kwadratów dopasuj dystrybuantę rozkładu cauchy'ego do dystrybuanty empirycznej.
Rozwiązanie
import scipy.optimize as so
def rho_cauchy(x,loc,scale):
return (np.pi*scale*(1.0+(x-loc)**2/(scale**2)))**(-1.0)
def F_cauchy(x,loc,scale):
return 0.5+np.arctan((x-loc)*1.0/scale)/np.pi
#losujemy 10000 liczb z rozkladu Caychyego o loc=1.23 i scale=2.0
x=2*np.random.standard_cauchy(10000)+1.23
N=len(x)
#METODA 1 - Dopasowanie metoda najmniejszych kwadratow do histogramu
#tworzymy histogram
hist,bins= np.histogram(x,bins=np.linspace(-20,20,61))
#dlugosc przedzialu histogramowania
przedzial=bins[1]-bins[0]
#normalizujemy histogram aby moc go porownac z gestoscia
hist=hist*1.0/len(x)/przedzial
#liczymy wsp. srodkow przedzialow histogramowania
xhist=bins[:-1]+0.5*przedzial
#definiujemy sume kwadratow
def squares((loc,scale)):
return sum([(rho_cauchy(xhist[i],loc,scale)-hist[i])**2 for i in range(len(hist))])
#szukamy minimum funkcja fmin
fit1=tuple(so.fmin(squares,(0,1)))
print 'wynik metody1 to '+str(fit1)
#ogladamy wynik
py.plot(xhist,hist)
xtest=np.linspace(-20,20,1001)
ytest=[rho_cauchy(a,*fit1) for a in xtest]
py.plot(xtest,ytest)
py.show()
#METODA 2 - Metoda najwiekszej wiarygodnosci
#definiujemy -funkcje wiarygodnosci
def L((loc,scale)):
return -sum([np.log(rho_cauchy(a,loc,scale)) for a in x])
#szukamy minimum
fit2=tuple(so.fmin(L,(0,1)))
print 'wynik metody2 to '+str(fit2)
#METODA 3 - dopasowanie dystrybuant
xx=sorted(x)
yy=np.linspace(0,1,N)
#definiujemy funkcje KS bedaca maksimum z roznicy miedzy dystrybuanta empiryczna a teoretyczna
def KS((loc,scale)):
return max([abs(F_cauchy(xx[i],loc,scale)-yy[i]) for i in xrange(N)])
#szukamy minimum
fit3=tuple(so.fmin(KS,(0,1)))
print 'wynik metody3 to '+str(fit3)
#ogladamy wynik
cut=100
py.plot(xx[cut:-cut],yy[cut:-cut])
xtest=np.linspace(-20,20,1001)
ytest=[F_cauchy(x,*fit3) for x in xtest]
py.plot(xtest,ytest)
py.show()
Zadanie - Data Container
Napisz klasę funkcja przyjmującą w konstruktorze parametry funkcji i posiadającąreprezentacje tekstową. Napisz dowolna funkcję dziedziczącą po klasie funkcja, której metoda call przyjmuje jeden argument i zwraca wartość funkcji dla podanego argumentu i parametrów podanych w konstruktorze. Napisz klasę DataContainer przyjmującą w konstruktorze dwie serie danych empirycznych o tej samej długości odpowiadające wspólrzędnym x i y. Obiekt klasy DataContainer powinien być wyposażony w metodę o nazwie fit, która przyjmuje jako argument obiekt klasy funkcja. Metoda fit powinna zwracać obiekt takiej samej klasy jak otrzymany w argumencie ale z parametrami dla których funkcja jest najlepiej (metodą najmniejszych kwadratów) dopasowana do przechowywanych w obiekcie danych.
Rozwiązanie
import numpy as np
import pylab as py
import time
import scipy.optimize as so
class funkcja(object):
def __init__(self,*args):
self.args=args
def __str__(self):
return 'to jest funkcja o nazwie '+self.__class__.__name__+' i argumentach '+str(self.args)
class liniowa(funkcja):
def __call__(self,x):
return self.args[0]*x*x+self.args[1]
class DataContainer(object):
def __init__(self,x,y):
self.x=np.array(x)
self.y=np.array(y)
self.n=len(x)
def __str__(self):
return '''to jest Data Container z danymi:
x[:10]:'''+str(self.x[:10])+'''
y[:10]:'''+str(self.y[:10])
def fit(self,funkcja):
parametry_poczatkowe=funkcja.args
def squares(parametry):
funkcja.__init__(*tuple(parametry))
return sum((map(funkcja,self.x)-self.y)**2)
parametry_dopasowane=so.fmin(squares,parametry_poczatkowe)
funkcja.__init__(*tuple(parametry_dopasowane))
return funkcja
#generujemy przykladowe xlist
xlist=np.arange(0,1,0.001)
#generujemy wartosci funkcji z szumem
f=liniowa(1.23,-0.73)
ylist=[f(x)+0.05*np.random.randn() for x in xlist]
d=DataContainer(xlist,ylist)
f=liniowa(1,2)
f=d.fit(f)
py.plot(d.x,d.y)
py.plot(d.x,map(f,d.x))
py.show()
