Curve fitting with NeuronC


Methods for curve fitting

To curve fit a model to a set of data, one runs the model with different sets of parameters to find the the output that best matches the data. For each parameter set, the model output is compared to the existing data with a comparison function, which is the difference between the model's output and the data. This function is minimized through a process of testing different parameter sets for the model. The difficulty in this process is finding an algorithm that finds the parameter set that produces the best match. It is difficult because in a nonlinear model the output (given a set of parameters) is not always easy to predict.

Nonlinear model fits

There are many methods for fitting parameter sets to a nonlinear model. One widely used method is the "Levenberg-Marquardt" method. It is useful when the model consists of a known continuous function that has computable first and second derivatives. This allows a very efficient "iterative" algorithm to find the best fit. The derivatives allow the iterative steps to be nearly optimal in size and direction.

The Levenberg-Marquardt least-squares fitting method is implemented in NeuronC. You can use the "lmfit()" procedure in either the interpreted or compiled versions, and a 2D compiled version is also available. A test program, "gaussfit.cc" is provided for the 2D version.

Another method is "downhill simplex" method. This is often useful when a simulation function can be evaluated but no first or second derivative is directly available. Since the simplex method does not rely on the derivatives of the function, it is not as computationally efficient, but it can work with a wider variety of nonlinear functions. Both of these methods work well only to find "local minima" and may get stuck in a minima which is not the global minimum for a function.

Random search methods

There are several methods for fitting nonlinear functions that use a random search method. These methods don't get stuck in local minima, and don't require knowing the derivatives of the simulated function, and in many cases can be faster than the simplex method. The reason is that an N-dimensional landscape in a typical neural circuit simulation can be very rugged, e.g. the voltage in a neuron is nearly linear with input for some combinations of parameters but for other combinations it fires action potentials which are very nonlinear. If the search includes a random component to it, the nonlinearities prevent the search from getting stuck in a local minimum and can find a global minimum faster than a purely directed search method.

Simulated annealing and stochastic search

One method is simulated annealing which finds a local minimum but has random noise added to the comparison function which allows it to jump out of a local minimum to be more efficient at finding a global minimum. One variation of this method adapts the "downhill simplex" method with different levels of noise depending on how close it gets to the minimum. Another method is stochastic search. This method is very general and does not use a directed search procedure. Instead it relies on random sampling to find the minima.

Both simulated annealing and stochastic search are implemented in NeuronC. To use these fitting methods:

  • 1) define a match function which computes the difference between simulation's output and a template of data values to match.

  • 2) define a simulation function called "runsim()" which runs the simulation and the match function and returns with the match function's value.

  • 3) Set up the test values and ranges for the free parameters and run the "ssrch()" search procedure which will run the search, calling the "runsim()" function to evaluate the match.

  • 4) For the implementation of simulated annealing look in "nc/tcomp/simann.n", and for stochastic search look in "nc/tcomp/stsrch.n". As an example of how to run these methods, see the scripts "nc/tcomp/testsa.n" and "nc/tcomp/rsrch".