\documentstyle[12pt]{article}
\baselineskip=18pt
\begin{document}
\title{\bf Fitting Curves by Simulated Annealing}
%
\author{T.J.P.Penna\\
Instituto de F\'{\i}sica, Universidade Federal
         Fluminense,\\
        Av. Litor\^anea, s/n\\ 24210-340 Niter\'oi, RJ,
         Brazil\\
         e-mail: gfitjpp@vmhpo.uff.br}
\date{\today}
\maketitle
\vfill
\begin{abstract}
Simulated annealing  is shown to be a powerful me\-thod to
fitting curves. This method is particularly appropriate to fit curves with
discontinous derivatives where the traditional methods fail. Recent
generalizations of the simulated annealing method enable easy and efficient
implementation on microcomputers.  
\end{abstract}
\vfill
\centerline{\em to appear in Computers in Physics }
\vfill
\newpage
Fitting observational data (by experiments or numerical simulations)
to a given  function is a  common task in science, sometimes to test
a theoretical model or as interpolating purposes, etc.  The most
popular method  consists  in minimizing the squared deviation
function
\begin{equation}
\varepsilon^2 = \sum_{i=1}^N \left\lbrack y_i - y(x_i;a_1\dots a_M)\right\rbrack^2 \label{chi}
\end{equation}
where $y_i$ $(i=1\dots N)$ are the data, $y(x)\equiv y(x,a_1 \dots
a_M)$ is the function to be fit, with the $M$ free parameters
$(a_1\dots a_M)$. Let us consider the data without error bars, for
the sake of simplicity. The errors $\sigma_i$ $ ( i=1\dots N )$ can
be introduced dividing each term on the right side of ( \ref{chi} )
by $\sigma_i$.  The specific methods to fitting curves frequently
demand a good estimate of the parameters to be fit\cite{nr}. In
addition, derivatives of the function have to been evaluated. This
fact restricts the application of those method to fit well-behaved
functions.  In this note, we show that a simulated annealing
approach\cite{kirk} can be useful in situations where the traditional
methods fails as, for example, for discontinous derivative functions.
Moreover, the implementation is easier and, due to recent
developments, the performance is also good.

The simulated annealing is widely used in the minimization of
nonconvex (multiple minima) functions -- as that presented in
eq.(\ref{chi}) -- as alternative to gradient descent methods. As far
as we know, it has not  been applied in fitting curves because other
methods are more efficient than the original version of the simulated
annealing.  Here we present some improvements in this approach that
enable it to work as an efficient and simple algorithm. The method is
implemented as follows: the variables are allowed to change (a jump
in the configuration space), according to a {\em visiting
distribution} that gives the probability of occurence of a jump of a
given size (we will return to this point in more detail later). If
this jump decreases the value of a cost function then the changes are
always accepted, just like gradient descent. However, if the jump is
in direction of increasing the cost function, the changes will be
accepted according to an  {\em acceptance probability}, depending on
an artificial {\em temperature}. This step avoids the trapping in
local minima. Then, the temperature is gradually cooled, according to a
recipe that depends on the choice of the visiting distribution, and
the process is repeat until the desired precision is reached.

Geman and Geman\cite{geman} showed that, in order to reach a global
minimum, if the visiting distribution is a {\em Gaussian} one, the
temperature must decrease logarithmically with time. This dependence
strongly discourages the use of simulated annealing as modeling data
algorithm.  Szu and Hartley\cite{szu} proposed to use a
Cauchy-Lorentz distribution (semi-local distribution) instead of the
Gaussian one (local distribution).  In this case, the temperature
decreases like the inverse of time, therefore the algorithm is more
efficient.  Recently, Tsallis and Stariolo\cite{gsa} proposed a new
visiting distribution that recovers both cited before as special
cases. A new  parameter $q_v \in \lbrack1,3\lbrack$ is introduced,
such $q_v=1$ recovers the classical simulated annealing (with
Gaussian distribution of jumps, sometimes called {\em Boltzmann
machine}), $q_v=2$ recovers the fast simulated annealing (Cauchy
machine).  The probability of trial value for a parameter $x$
changing from $x_t$ to $x_{t+1}$ is\cite{gsa}
\begin{equation}
g_{q_V}(\Delta x_t) \propto \frac{T(t)^{-D/(3-q_v)}} {
\left\{1+(q_v-1)(\Delta x_t)^2/ T(t)^{2/(3-q_v)}\right\}
^{[1/(q_v-1)] + [(D-1)/2]} }
\end{equation} 
where $T(t)$ is the artificial temperature at time $t$, $D$ is the
dimension of the function to be minimized and $\Delta x_t= x_{t+1} -
x_t$ is the size of the jump. For arbitrary $D$, the temperature
should be cooled as
\begin{equation}
T(t) = T(1) \frac {2^{q_v-1}-1} {(1+t)^{q_v-1}-1}.\label{cool}
\end{equation}  
Therefore the temperature decreases as a $(q_v-1)$-power of time,
consequently in a much more efficient way than the earlier
strategies.

Now, let us focus the discussion on the acceptance probability as
function of temperature. The most frequent used is the so-called
Metropolis algorithm\cite{metro}. Since the temperature is an
artificial parameter the only condition to be respected is  whether
the acceptance probability leads to correct equilibrium states.
Within the framework of a generalized statistical mechanics
introduced by Tsallis\cite{tsallis}, we can define a generalized
version of the Metropolis algorithm\cite{tsp}: the new value of the
function to be minimized is accepted with a probability given by
\begin{equation}
p=\mbox{min}\left(1, ( 1 - (1-q_a)\left( y(x_{t+1}) - y(x_t)
\right)/T(t) )^{(1/  (1-q_a)}\right) \label{gma}
\end{equation}
where $q_a$ is a new parameter.  We recover the Metropolis algorithm
in the $q_a\rightarrow 1$ limit. The temperature for the acceptance
probability can be different from the temperature for the visiting
distribution. Here, we made them equal, for simplicity.
Eq.(\ref{gma}) also leads to equilibrium states for all values of
$q_a$\cite{gsa,tsp}.  The most important feature of the
eq.(\ref{gma}) -- for our purposes -- is that approximate solutions
are found even at higher temperatures (for negative values of $q_a$)
than had we used the Metropolis algorithm. Therefore, starting  from
a high temperature we {\em do not need} good initial values.

In fig.1, we present a complete code written in C where we try to fit
a Gaussian function centered at origin. The variance and the
amplitude are the parameters to be fitted. In order  to show this
method is useful to fitting curves, we present a situation where
there are local minima. In the fourth line of the code, we define the
function to be fitted (it will used also to generating the input
data): two combinated Gaussian functions, with dispersion $A$ for
negative values of $X$, and dispersion $C$ for the positive ones. The
amplitude is given by parameter $B$. For  input data we have chosen
$A=2.0$, $B=1.0$ and $C=2.15$. In this test case, we will try to fit
a single Gaussian function (frequently in modeling data we do not
know, {\it a priori}, the actual distribution of data, so this
example should not be considered as a biased one).  The error
landscape is strongly dependent of the distribution of the input
data, i.e., whether we have more data either in the positive or in the
negative section of axis $X$.  In this test, the input data were
equally spaced in the $X$-axis in the interval $\lbrack
-1.7,2.3\rbrack$.  In fig.2, one can see how the variance evolves in
time when we started from a low temperature (30 times smaller than
initial value of the error). The system is initially trapped in a local minimum ($A=2.0$, in this case) but, after about 100 iterations, it goes to the global minimum ($A=
2.1137\pm 0.0005$). In this we can also see how accurate is the value of $A$,
checking the fluctuactions.  In this example the existence of
multiple minima turns evident but the temporary trapping  can be
avoided starting from higher temperatures. In this test, we follow a different approach from the one suggested by Stariolo and Tsallis\cite{gsa}:  we have allowed to change only one  parameter by step. In our computer tests, this approach has been seen as more efficient than the original one in finding the global minimum.
It is straightforward to extend this version for more than two parameters, just changing the value of $D$. Finally, we show the best fit for the test case, in fig.3.
 
In summary, we have applied successfully the generalized simulated
annealing to fitting curves. The fast cooling scheme and a
generalization of the Metropolis algorithm make this method
competitive with the traditional others, with additional advantages
as it avoids derivatives evaluation and good initial values. This
method was used\cite{prl} to fit more complex function as a
L\'evy-distribution
\begin{equation}
P(I,\psi,\gamma) = {1\over \pi} \int_0^\infty \exp(-\gamma q^\psi)
\cos(qI) dq
\label{plevy}
\end{equation}
to experimental data from leaky faucets, with $\gamma$ and $\psi$ as
free parameters and the integral was evaluated numerically\cite{nr}.

\noindent{\bf  Acknowledgments}

Financial support from the Brazilian Agencies CNPq and FINEP is
gratefully acknowledged.

\newpage
\begin{thebibliography}{99}
\bibitem{nr} W.H. Press, S.A. Teukolsky, W.T. Vetterling and B.P.
Flannery, {\it Numerical Recipes in C}, Cambridge University Press,
NY (1992).

\bibitem{kirk} S. Kirkpatrick, C.D. Gellat and M.P. Vecchi, {\it Science} {\bf 220}, 671 (1983).\\
S. Kirkpatrick, {\it J.Stat.Phys} {\bf 34}, 975 (1984).

\bibitem{geman} S.Geman and D. Geman, IEEE Trans. Patt. Anan. Mach. Int., PAMI-6 {\bf 6}, 721 (1984).

\bibitem{szu} H.Szu and R. Hartley, {\it Phys.Lett.} {\bf A122}, 157 (1987).

\bibitem{gsa}
D.A. Stariolo and C.Tsallis, in {\em Annual Review of Computational
Physics II}, ed. D.Stauffer, World Scientific Po., Singapore (1994).

\bibitem{metro} N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller and E. Teller, {\it J.Chem.Phys.} {\bf 21}, 1087 (1953).

\bibitem{tsallis} C. Tsallis, {\it J.Stat.Phys. } {\bf 52}, 479 (1988).

\bibitem{tsp} T.J.P. Penna,  {\it Phys.Rev.} {\bf E} (1994).

\bibitem{prl} T.J.P.Penna, P.M.C. de Oliveira, J.C. Sartorelli, W.M.Gon\c{c}alves and R.D. Pinto, {\em ``Long Range Anticorrelations: Is the Healthy Heart a Leaky Faucet ?'' }, preprint (1994).

\end{thebibliography}

\newpage
\noindent{\Large\bf Figure Captions}

{\bf Fig.1} A code written in C, for fitting  a gaussian function.

{\bf Fig.2} Time evolution of the variance.

{\bf Fig.3} Input data (squares) fitted by a Gaussian function (solid line). 
\end{document}
