Nonlinear function fitting in Matlab using Fminuit

Japanese version is here.

Please follow this instruction to get my software on your computer.

A good ending of data analysis is making a model that explains your data. Typically it requires nonlinear function fitting. In Matlab, unless you have toolboxes, you will be left with ‘fminsearch’ function with limited capability. It is nice if you can do more sophisticated nonlinear fitting without a toolbox. In this post, I’ll introduce you to a very nice external function, ‘fminuit’, that does exactly that.

Fminuit is based on CERN‘s function fitting library Minuit. You can see a full explanation of what Fminuit is on its website. Briefly, Minuit is a program that is used in CERN for decades in order to fit high-energy physics models to petabytes of data collected there. It is very robust, efficient, feature-rich, and highly controllable. Fminuit is a Matlab wrapper of Minuit. My git repository includes only a binary version of fminuit for 64-bit Intel Macs. Windows or Linux binary are available from the Fminuit page, so please download and put them to your Matlab path. (You can also compile it.)

Fminuit itself is already pretty strong software, but I further wrapped it by using my ErrorNum class to make it even more user-friendly. Let’s take a look at the beginning of the help of my ‘MinuitFitEN’ function.

  function res = MinuitFitEN(func, initp, x, y, varargin)
 
     func - fit funciton handle func(p, x)
     initp - (nPar vec) estimated inital parameters
     x - (any) 2nd argument of func
     y - (any) ErrorNum data to fit including error
 
     -- from here optional parameters (default)
     'Extradim' - (scalar) extra degree of freedom to subtract (0)
     'LowerLimit' - (nPar vec) lower bound of the parameters (false)
     'UpperLimit' - (nPar vec) upper bound of the parameters (false)
     'StepSize' - (nPar vec) step sizes of the parameters (false)
     'minfunc' - (String) function to minimize ('PoissonNLL')
     'command' - (String) command to give to fminuit
                 ('set pri -1; scan; minimize; simplex; improve; minos')

You can run a nonlinear function fit by writing as simple as the following example. First, let’s construct the data to fit.

x = 0:0.2:3;
n_obs = 30; % number of observation
sigma = 5; % STD of the data

% Through the fit, you need to estimate these variables
base = 2; % base of the power function
offset = 3; % offset from 0

% Construct your data based on the parameters above
% Note the last term is adding errors in observations.
y_obs = repmat(base.^x, n_obs, 1) + offset + sigma * randn(n_obs, length(x));
y_err = ErrorNum.create(y_obs, 1); % ErrorNum.create calcualtes SEM

figure(101);
plotEN(x, y_err, '.', 'capsize', 0); % EN functions are for the ErrorNum class.

This yields a figure like this. As expected, you get data points that roughly obey the power function, with some fluctuations. Now, let’s fit a function to it.

powfunc = @(p, x) p(1).^x + p(2); % this is your model function to fit the data.

% you need to guess the initial parameters
base_guess = 3;
offset_guess = 1;
initial_p = [base_guess, offset_guess];
fitresult = MinuitFitEN(powfunc, initial_p, x, y_err, 'minfunc', 'Chi2');

In the first line, I’m defining an anonymous function as my model. If you have a .m file, you can get the handle by doing something like ‘@function_name’. Now fitting is done. The last two arguments for MinuitFitEN indicate the type of fitting (or a function to minimize). This time, I am doing a chi-square fit, so I’m choosing ‘Chi2’. As I’m doing a lot of analysis on spikes of neurons, the default is set to ‘PoissonNLL’, or Negative Log-Likelihood function for the Poisson distribution. In applications where ANOVA works (like this example), ‘Chi2’ should be the choice.

Let’s take a look at the result of parameter estimation.

>> fitresult.p
ans = 
ErrorNum:
2.07 ± 0.06,  2.8 ± 0.3

The output of the MinuitFitEN function is a structure that contains the estimated parameters with errors. The errors are one sigma of ANOVA, so you have a ~68% chance of having your true value within this error. The first value should be 2, but it’s slightly off. The second value is nicely within our expectation, 3. The fit seems to be working! It is very nice to always have the error attached to your results. The output is again, an instance of the ErrorNum class, so we get a reasonable text output without adjusting anything.

MinuitFitEN also calculates goodness-of-fit statistics. They are stored in these variables.

>> fitresult.redchi2
ans =
    1.5140

>> fitresult.redchi2sig
ans =
    0.9033

‘redchi2’ stands for ‘reduced chi-squared’, I wouldn’t give a detailed explanation of this value. Please refer to Wikipedia for its definition and usage. Here, I just say that if this value is close to 1, you are doing a reasonable fitting. If it is too small, you are overfitting (or overestimating the error); if it is too large, you are underfitting (or underestimating the error).

‘redchi2sig’ gives you the significance of the goodness of your fit. It is like a p-value for the goodness-of-fit. For example, if it is larger than 0.99, you have a reduced chi-squared value that a perfect fit would encounter only 1% of the chance. To say that your model explained the data well, you want to have a reasonable value (e.g. 0.05 < redchi2sig < 0.95).

These are important statistics, and I extensively used these values in the last two papers (Ito et al., 2017; Ito et al., 2020) that I have written. (These papers have real-world examples of function fitting.) However, these are also difficult concepts to grasp if you are a beginner of function fitting. In future posts, I will demonstrate more cases to show how to use these values. I will also demonstrate more advanced usage of this function (such as setting parameter boundaries, using more complicated functions, etc…), but for now, please refer to its help to see how to use it.