Függvényillesztés Monte Carlo módszerrel mérésleírás

Az alábbi leírás pdf formátumban is letölthető. Példa adatsor: Data.txt.

Bevezetés

A Monte Carlo algoritmusok szerteágazó alkalmazási területei közül a függvényillesztés a számítógépes gyakorlat témája. A leírás áttekinti a szükséges elméleti hátteret, valamint segítséget ad a konkrét megvalósításhoz.

Elméleti háttér

Legyen az illesztendő adatsor $\setbox0\hbox{$ Y_\textrm{orig}(X_\textrm{i})$}% \message{//depth:\the\dp0//}% \box0%$, ahol $\setbox0\hbox{$i=1\ldots N$}% \message{//depth:\the\dp0//}% \box0%$. Az illesztett adatsor hasonlóan $\setbox0\hbox{$Y_\textrm{fit}(X_\textrm{i}, \vec{A})$}% \message{//depth:\the\dp0//}% \box0%$, ahol $\setbox0\hbox{$\vec{A}=(A_1,A_2 \ldots A_\textrm{M})$}% \message{//depth:\the\dp0//}% \box0%$ a paraméterek listája, amelyek hangolásával érjük el, hogy az illesztett adatsor rásimuljon az eredetire. Ez a következő energiakifejezés minimalizálását jelenti:

$\[E(\vec{A})=\sum_{i=1}^{N} \left| Y_\textrm{orig}(X_\textrm{i}) - Y_\textrm{fit}(X_\textrm{i}, \vec{A})\right| \label{MCEnergyEq} \]$

A minimalizálást Monte Carlo--Metropolis algoritmussal végezzük: Változtassuk meg a paramétereket véletlenszerűen: $\setbox0\hbox{$\vec{A}' = \vec{A}+ \vec{\delta A}$}% \message{//depth:\the\dp0//}% \box0%$. Ezzel a fentiek szerint kiszámolt energia is változik: $\setbox0\hbox{$E(\vec{A}) \rightarrow E(\vec{A}')$}% \message{//depth:\the\dp0//}% \box0%$. Az új paramétereket elfogadjuk, ha közelebb kerültünk a minimumhoz, tehát $\setbox0\hbox{$E(\vec{A}') < E(\vec{A})$}% \message{//depth:\the\dp0//}% \box0%$. Célunk elkerülni, hogy az algoritmus befagyjon egy lokális minimumba, ezért $\setbox0\hbox{$p=\textrm{exp}(-\beta(E(\vec{A}')-E(\vec{A})))$}% \message{//depth:\the\dp0//}% \box0%$ valószínűséggel dolgozunk tovább az új paraméterekkel, ha $\setbox0\hbox{$E(\vec{A}') > E(\vec{A})$}% \message{//depth:\the\dp0//}% \box0%$. Vegyük észre, hogy $\setbox0\hbox{$\beta \rightarrow \infty$}% \message{//depth:\the\dp0//}% \box0%$ esetben csak lefelé léphetünk, míg a $\setbox0\hbox{$\beta \rightarrow 0$}% \message{//depth:\the\dp0//}% \box0%$ határátmenetben véletlenszerű mozgást kapunk a paramétertérben. Ez adja a fizikai jelentését az algoritmusban szereplő $\setbox0\hbox{$\beta$}% \message{//depth:\the\dp0//}% \box0%$ értéknek: megfeleltethető egy, a paramétertérben mozgó részecske $\setbox0\hbox{$(kT)^{-1}$}% \message{//depth:\the\dp0//}% \box0%$ inverz hőmérsékletének.

Implementáció

A fentiek megvalósításához a következő felépítést érdemes követni:

1. Készítsünk globális tömböt az eredeti, valamint az illesztett adatsornak;
2. Hasonlóan globális változókban tároljuk az inverz hőmérsékletet, az illesztési paramétereket, valamint azok lépésközét;
3. Deklaráljunk két véletlenszám-generátort: egyet a paraméterek megváltoztatásához és egyet az $\setbox0\hbox{$\exp(-\beta \Delta E)$}% \message{//depth:\the\dp0//}% \box0%$ valószínűségű lépéshez;
4. Hozzunk létre globális változókat a megváltoztatott paramétereknek, valamint egy tömböt a megváltoztatott paraméterekkel generált illesztőfüggvénynek.

Az illesztés inicializálásához hozzuk létre a szükséges változókat, és olvassuk be a kezdőértékeket a TextBoxokból.

Az illesztést az alábbiak szerint végezzük:

1. Számítsuk ki a $\setbox0\hbox{$E(\vec{A})$}% \message{//depth:\the\dp0//}% \box0%$ szerinti energiát az aktuális változókkal;
2. Változtassuk meg az egyik paramétert a megadott határok között egyenletes eloszlással;
3. Számítsuk ki az így megváltozott energiát;
4. A régi és új energia különbsége alapján döntsünk arról, hogy elfogadjuk-e a változást: ha az új érték kisebb, akkor minden esetben, ha nem, akkor sorsoljunk egy véletlenszámot, és $\setbox0\hbox{$\exp(-\beta \Delta E)$}% \message{//depth:\the\dp0//}% \box0%$ valószínűséggel mégis elfogadjuk a változást;
5. Ha szükséges, írjuk ki az új értéket a TextBoxba, és ábrázoljuk a megváltozott illesztőfüggvényt;
6. A fentiek szerint járjunk el a többi paraméter esetében.

Feladatok

Készítsünk grafikus függvényillesztő felületet a data.dat file-ban található 500 elemű adatsor illesztéséhez a

$\[ y(t)=Ae^{-t/\tau}\sin(2\pi t/T)\]$

függvénnyel, ahol $\setbox0\hbox{$A,\tau,T$}% \message{//depth:\the\dp0//}% \box0%$ a keresett paraméterek. A cél eléréséhez az alábbiak szerint haladjunk:

1. Olvassuk be és ábrázoljuk a data.dat állomány tartalmát!
2. Hozzunk létre TextBoxokat, amelyekkel állíthatóak a függvényillesztés kiinduló paraméterei, azok lépésköze, valamint az inverz hőmérséklet!
3. Hozzunk létre gombokat, amelyekkel elindíthatjuk a függvényillesztést: legyen lehetséges újraindítani az illesztést a TextBoxokban megadott kezdeti paraméterekkel, valamint 1, illetve 100 illesztési ciklust elindítani!
4. A program jelenítse meg az aktuális paramétereket, az $\setbox0\hbox{$E(\vec{A})$}% \message{//depth:\the\dp0//}% \box0%$ szerinti energiát, valamint az illesztett függvényt!

A jegyzőkönyvben szerepeljen a program felépítésének leírása, valamint a probléma diszkussziója. Vizsgáljuk meg a kezdeti paraméterek, az inverz hőmérséklet, valamint a paraméterek lépésközeinek hatását az illesztés gyorsaságára és jóságára, azaz adjuk meg, hogy átlagosan hány lépés után illeszkedik a generált függvény az eredetire, valamint mekkora a jellemző minimális energia! Opcionálisan próbáljuk ki, hogy lehetséges-e jobb eredményt elérni, ha $\setbox0\hbox{$\beta$}% \message{//depth:\the\dp0//}% \box0%$-t, illetve a lépésközöket dinamikusan változtatjuk! A dokumentációhoz csatoljuk a program forráskódját!

Syllabus in English

The exercise sample data is available here: Data.txt.

Introcution

Monte Carlo algorithms have a widespread use in many fields, this exercise focuses on fitting a function to a series of data (such as what would be obtained by an experimental measurement). The exercise description presents the necessary theoretical background, and provides indications for the concrete implementation.

Theoretical background

Let the original series of points we desire to fit, be denoted as $\setbox0\hbox{$ Y_\textrm{orig}(X_\textrm{i})$}% \message{//depth:\the\dp0//}% \box0%$, where $\setbox0\hbox{$i=1\ldots N$}% \message{//depth:\the\dp0//}% \box0%$. The fitted data series is $\setbox0\hbox{$Y_\textrm{fit}(X_\textrm{i}, \vec{A})$}% \message{//depth:\the\dp0//}% \box0%$, where $\setbox0\hbox{$\vec{A}=(A_1,A_2 \ldots A_\textrm{M})$}% \message{//depth:\the\dp0//}% \box0%$ is the list of parameters that are varied to obtain the best correlation between the fitted series and the original series. The best fit is to be regarded as a minimization of the following "energy" expression:

$\[E(\vec{A})=\sum_{i=1}^{N} \left| Y_\textrm{orig}(X_\textrm{i}) - Y_\textrm{fit}(X_\textrm{i}, \vec{A})\right| \label{MCEnergyEq} \]$

The minimization of energy is obtained through the Monte Carlo - Metropolis algorithm: First, alter the fitting parameters with a random factor: $\setbox0\hbox{$\vec{A}' = \vec{A}+ \vec{\delta A}$}% \message{//depth:\the\dp0//}% \box0%$. Using the new randomly varied parameters, a new energy value is obtained: $\setbox0\hbox{$E(\vec{A}) \rightarrow E(\vec{A}')$}% \message{//depth:\the\dp0//}% \box0%$. The new parameters should be accepted and kept, if the new energy is closer to the energy minimum, that is, $\setbox0\hbox{$E(\vec{A}') < E(\vec{A})$}% \message{//depth:\the\dp0//}% \box0%$. The key to the algorithm is to avoid becoming stuck in a local minimum, thus we use an "annealing" method, where with $\setbox0\hbox{$p=\textrm{exp}(-\beta(E(\vec{A}')-E(\vec{A})))$}% \message{//depth:\the\dp0//}% \box0%$ probability we accept the new parameters, even if $\setbox0\hbox{$E(\vec{A}') > E(\vec{A})$}% \message{//depth:\the\dp0//}% \box0%$. Observe that using $\setbox0\hbox{$\beta \rightarrow \infty$}% \message{//depth:\the\dp0//}% \box0%$ the algorithm can only step downwards in energy, while in the limit $\setbox0\hbox{$\beta \rightarrow 0$}% \message{//depth:\the\dp0//}% \box0%$ the walk in the parameter space is entirely random. This gives a physical meaning to the $\setbox0\hbox{$\beta$}% \message{//depth:\the\dp0//}% \box0%$ value: it can be regarded as particle moving in the parameter space with $\setbox0\hbox{$(kT)^{-1}$}% \message{//depth:\the\dp0//}% \box0%$ inverse temperature.