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

A Fizipedia wikiből

See the English version at the bottom of the page.

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


Tartalomjegyzék

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 N (TextBoxban megadható) 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!
  5. Célszerű egy külön grafikonon ábrázolni az \setbox0\hbox{$E(\vec{A})$}% \message{//depth:\the\dp0//}% \box0% előző értékeit lépésszám függvényében. Ezen szemmel könnyen áttekinthető az algoritmus konvergenciája.

A jegyzőkönyvben szerepeljen a program felépítésének leírása, valamint a probléma diszkussziója. Vizualizáljuk az \setbox0\hbox{$E(\vec{A})$}% \message{//depth:\the\dp0//}% \box0% skalárteret alacsonyabb dimenziós (1D és 2D) metszeteken keresztül. 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. Grafikonon ábrázoljuk, az inverz hőmérséklet néhány nagyságrendnyi függvényében, 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 (átlagoljuk ki több illesztési próbálkozásból). 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.

Implementation


To create the program with the above considerations, follow these steps:

  1. Create a global variable array for both the original and fitted data series.
  2. Also in global variables, store the inverse temperature, the fitting parameters, and their maximum random increment;
  3. Create two random number generators: one for the parameter variation and one for the \setbox0\hbox{$\exp(-\beta \Delta E)$}% \message{//depth:\the\dp0//}% \box0% probability acceptance step;
  4. Create global variables to temporarily store the randomly altered parameter set and a data series generated by this set. This is needed so that reverting to the previous parameter set and series is possible in case a new parameter set is not accepted.

To initiate the fitting, initialize the required variables, and read the starting values from TextBoxes.

The steps in the fitting algorithm are the following:

  1. Calculate \setbox0\hbox{$E(\vec{A})$}% \message{//depth:\the\dp0//}% \box0% energy using the current fitted series;
  2. Alter one of the parameters within the specified limits with a uniform random distribution;
  3. Calculate the new energy;
  4. Based on the difference between the new and previous energy, decide to accept the parameter change: If the new value is smaller, then accept it; if not, then get a new random value between 0-1, and with \setbox0\hbox{$\exp(-\beta \Delta E)$}% \message{//depth:\the\dp0//}% \box0% probability accept it anyway; otherwise discard it;
  5. If necessary, display the new values in a TextBox and display the altered fitting function on graph;
  6. Repeat these steps for every other parameter as well.


Tasks

Create a graphical interface for controlling the fitting program. Use the provided data.txt as an input, having 500 elements in the data series. The fitting is to be attempted using the
\[ y(t)=Ae^{-t/\tau}\sin(2\pi t/T)\]
function where \setbox0\hbox{$A,\tau,T$}% \message{//depth:\the\dp0//}% \box0% are the parameters to find. To obtain the best fit, proceed the following way:
  1. Read the data.txt file contents into your program's memory and display its contents!
  2. Create TextBoxes, for the user to set the starting parameters of the function, their random variation interval, and the inverse temperature!
  3. Create buttons to start the fitting algorithm: let it be possible to restart the fitting using newly introduced TextBox values, to execute 1 step of the algorithm, or to execute a user-defined number of steps (given in a separate TextBox).
  4. The program should display the current actual parameters, the energy \setbox0\hbox{$E(\vec{A})$}% \message{//depth:\the\dp0//}% \box0%, and on graph the fitted series in a separate color than the original data series!
  5. It is advisable to display, in a separate graph, the history of \setbox0\hbox{$E(\vec{A})$}% \message{//depth:\the\dp0//}% \box0% energy as a function of fitting steps done by the algorithm. This makes it possible to tell when the algorithm has converged to a minimum.

In the report, describe the workings of your program, and discuss the problem. Visualize the \setbox0\hbox{$E(\vec{A})$}% \message{//depth:\the\dp0//}% \box0% scalar field through lower-dimensional (1D and 2D) cross-sections. Describe, with examples given, how your program reacts when different starting parameters are given, when various inverse temperatures are used, or how the random variation interval of parameters affects the output of fitting. You should describe the performance of your program by considering the number of steps necessary to find the minimum with certain starting conditions, as well as compare the minimum energies obtained, to be able to tell if some conditions prevented your algorithm from finding the global minimum energy. Optionally, try to check if you can get a better result (faster fit) if \setbox0\hbox{$\beta$}% \message{//depth:\the\dp0//}% \box0% and and variation intervals are dynamically tuned! Attach the source code to your documentation!