How to model a time series data with a custom function?

I'm trying to fit an ex-gaussian function to a time series of data. I tried to use the fit function of MATLAB. I don't know what I'm doing wrong, I'm not finding the best fitting parameters of my function.

Here is my code and an image of the results (plot on the right).

% Ex-Gaussian model
  exGaussian = @(mu,sigma,tau,a0,a1,time) a0 + a1/2 .* exp((1./(2*tau)) .* (2.*mu+(sigma.^2./tau)-(2.*time))).*...
     erfc((mu+(sigma.^2./tau)-time)./sqrt(2*sigma));

Fs = 20; %sampling frecuency
t = 0:1/Fs:5; %timestamp in seconds

%Parameters for the example data to fit
mu = 2; 
sigma = .1; 
tau = 1;
a0 = 1;
a1 = 19;

%Example data
y  = exGaussian(mu,sigma,tau,a0,a1,t);

figure(1),clf
subplot(1,2,1)
plot(t, y), grid off
title('Example of data')

% Select the options to fit the model
opts = fitoptions('Method','NonlinearLeastSquares',...
     'Lower',[0 0 0 5 10],'Upper',[10 2 10 30 50],...
     'Startpoint',[1 0.1 0.1 1 1]);
 
exGfit = fittype('exGaussian(mu,sigma,tau,a0,a1,t)',...
    'coefficients',{'mu','sigma','tau','a0','a1'},...
    'independent','t','options',opts);

subplot(1,2,2)
[fitResult,  fitStats, fitInfo] = fit(t',y',exGfit);
plot(fitResult,'-r',t,y,'.-b');

[Walter Roberson]

Guassians are difficult to fit -- some of the parameters need to be pretty close to the correct value or else you are very likely to get results that are no-where near correct.

That said, in this case it doesn't even do a great job when given the actual parameters as the starting point :(

% Ex-Gaussian model
  exGaussian = @(mu,sigma,tau,a0,a1,time) a0 + a1/2 .* exp((1./(2*tau)) .* (2.*mu+(sigma.^2./tau)-(2.*time))).*...
     erfc((mu+(sigma.^2./tau)-time)./sqrt(2*sigma));

Fs = 20; %sampling frecuency
t = 0:1/Fs:5; %timestamp in seconds

%Parameters for the example data to fit
mu = 2; 
sigma = .1; 
tau = 1;
a0 = 1;
a1 = 19;

%Example data
y  = exGaussian(mu,sigma,tau,a0,a1,t);

figure(1),clf
subplot(1,2,1)
plot(t, y), grid off
title('Example of data')

% Select the options to fit the model
opts = fitoptions('Method','NonlinearLeastSquares',...
     'Lower',[0 0 0 5 10],'Upper',[10 2 10 30 50],...
     'Startpoint',[mu, sigma, tau, a0, a1]);
 
exGfit = fittype(exGaussian,...
    'coefficients',{'mu','sigma','tau','a0','a1'},...
    'independent','time','options',opts);

subplot(1,2,2)
[fitResult,  fitStats, fitInfo] = fit(t',y',exGfit)

fitResult = 
     General model:
     fitResult(time) = a0+a1/2.*exp((1./(2*tau)).*(2.*mu+(sigma.^2./tau)-(2.*time)
                    )).*erfc((mu+(sigma.^2./tau)-time)./sqrt(2*sigma))
     Coefficients (with 95% confidence bounds):
       mu =       2.358  (1.745, 2.97)
       sigma =     0.07661  (-0.009826, 0.163)
       tau =      0.3223  (-0.0268, 0.6714)
       a0 =           5  (fixed at bound)
       a1 =       16.62  (13.08, 20.17)

fitStats = struct with fields:
           sse: 607.8905
       rsquare: 0.5742
           dfe: 97
    adjrsquare: 0.5611
          rmse: 2.5034

fitInfo = struct with fields:
           numobs: 101
         numparam: 5
        residuals: [101×1 double]
         Jacobian: [101×5 double]
         exitflag: 3
    firstorderopt: 1.4248
       iterations: 12
        funcCount: 78
     cgiterations: 0
        algorithm: 'trust-region-reflective'
         stepsize: 0.0030
          message: 'Success, but fitting stopped because change in residuals less than tolerance (TolFun).'

plot(fitResult,'-r',t,y,'.-b');

Answer 1

[Mathieu NOE]

hello

a quick and dirty trial with the poor's man solution ( with fminsearch)

the only trick to avoid an error (because with negative sigma, sqrt(sigma) is complex), I modified the function :

sqrt(2*abs(sigma)))

not 100% scientific maybe regarding optimizer usage, but it's seems to work not too bad

made a few trials with different IC, not all will give a good fit but with a liitle effort we can probably refine the IC range

% Ex-Gaussian model
  exGaussian = @(mu,sigma,tau,a0,a1,time) a0 + a1/2 .* exp((1./(2*tau)) .* (2.*mu+(sigma.^2./tau)-(2.*time))).*...
     erfc((mu+(sigma.^2./tau)-time)./sqrt(2*abs(sigma)));


Fs = 20; %sampling frecuency
t = 0:1/Fs:5; %timestamp in seconds

%Parameters for the example data to fit
mu = 2; 
sigma = .1; 
tau = 1;
a0 = 1;
a1 = 19;

%Example data
y  = exGaussian(mu,sigma,tau,a0,a1,t);
figure(1),clf
plot(t, y), grid off
title('Example of data')


% curve fit using fminsearch
obj_fun = @(p) norm(exGaussian(p(1),p(2),p(3),p(4),p(5),t)-y);
p0 = [0.5 0.01 0.5 y(1) max(y)]; % IC guess for mu,sigma,tau,a0,a1

sol = fminsearch(obj_fun, p0);
yfit = exGaussian(sol(1),sol(2),sol(3),sol(4),sol(5), t);
plot(t,yfit,'-',t,y,'r .','MarkerSize', 20);

title('Data', 'FontSize', 20)
ylabel('Amplitude', 'FontSize', 20)
xlabel('t', 'FontSize', 20)
legend('curve fit','raw data');

[María ML]

Thanks! That simple change in the equation will help me a lot.

[Mathieu NOE]

My pleasure !

[George]

Thanks so much for this Mathieu, it works for me too. Pleased to have found it.

I'm trying to understand how it works as well. The Ex-Gaussian formula is a bit different from one I have seen before. In particular I am not sure what the parameters a0 and a1 do, so if you had some more information on that I would be really grateful.

Many thanks,

George

[Mathieu NOE]

hello @George

a quick search on internet , found this wiki publication

Exponentially modified Gaussian distribution - Wikipedia

so , I believe a0 is simply a dc constant that has been added to cope with curves that start at y different from zero (vertical shift) , and a1 is same as lambda

Answer 2

[Prateek]

Hi Ana,

A good way to fit ex-Gaussian function is to use the "interpolant" option in the curve fitting toolbox. I was able to obtain a good fit for the data providede by you:

Hope this helps.

Prateek