Problem in data fitting using nonlinear regression model

Hi,

I am having issues to fit data with the nonlinear regression model using 'fitnlm' command. Here attached my code and the model data.

Following errors appear-

Error using nlinfit>checkFunVals (line 649)

The function you provided as the MODELFUN input has returned Inf or NaN values.

Any suggestions?

%Normalized transmission mode pump-probe data fitting using Hydrodynamic
%model
close all
clear all
clc
%Matlab code written by Subhajit Bej and edited Sudip Gurung
%Dated: 04/10/2020
%Edited by Subhajit Bej on 18.04.2020
%Edited and corrected by Subhajit Bej on 22.04.2020
%Edited by Subhajit Bej on 29.04.2020

% AA=readmatrix('ITO_power_K(x), L(10GW), M(21GW).xls');%read data from excel (97-2003) workbook and save in matrix format
% x=AA(:,11); %import your xdata. Contains centered position pump-probe delay (tau) data. Eleventh column in the excel sheet
% y=AA(:,15); %import your ydata. Contains normalized transmission data. Maximum/ minimum position must match with tau=0. Fiftennth column in the excel sheet
% x=xlsread('ITO_power_K(x), L(10GW), M(21GW).xls','K4:K202');
% y=xlsread('ITO_power_K(x), L(10GW), M(21GW).xls','O4:O202');

% x=xlsread('ITO_power_K(x), L(10GW), M(21GW)_test.xls','K4:K130');
y=1+xlsread('ITO_power_K(x), L(10GW), M(21GW)_test.xls','T4:T130');
y=y';
A=627;
N=length(y);%Number of points
x=linspace(-A,A,N);%alternatively this command can also be used to define 't' vector
% x=x';

%Finding the x position of the peak for the experimental data
xIndex0=find(y==max(y));
x0=x(xIndex0);

%Setting the optimization type and parameters
opts=statset('fitnlm');
opts.RobustWgtFun ='cauchy';% A weight function for robust fitting. Valid functions
%                        are 'bisquare', 'andrews', 'cauchy', 'fair', 'huber',
%                        'logistic', 'talwar', or 'welsch'. Default is '' (no
%                        robust fitting). Can also be a function handle that
%                        accepts a normalized residual as input and returns
%                        the robust weights as output.
opts.Robust='on';
opts.TolFun=1e-25;
opts.MaxIter=10000;
opts.MaxFunEvals=5000;
opts.Tolx=1e-25;
opts.Display='iter';

%Fitting the experimental data using Nonlinear regression technique.
%beta0=[50;90;6;1;1];%Initial guess of the parameters to be fitted. It is a good idea to define the parameters here.
beta0=[24;50;120;4];%Initial guess of the parameters to be fitted. It is a good idea to define the parameters here.
% beta(1)= probe pulse duration;
% beta(5)= amplitude of probe pulse;
% beta(4)= amplitude of response time;
% beta(3)= thermalization time;
% beta(2)= recombination time;

mdl=fitnlm(x,y,@hydrodynamic,beta0,'Options',opts);
beta=mdl.Coefficients.Estimate;
r2=mdl.Residuals.Raw;

ydash=hydrodynamic(beta,x);

%plots
plot(x,y,'k--o','LineWidth',1);
hold on;
plot(x,ydash,'r-','LineWidth',2.4);
hold off;
xlabel('\tau [fs]');
ylabel('\Delta T/\Delta T_{max}');
legend('experimental data','fitted plot');


function G=hydrodynamic(beta,x) %For pump-probe data fitting.
%Probe pulse
% Probe=gauspuls(t,fc,tauprobe);%Probe pulse with temporal Gaussian profile
% Probe=gauspuls(t,fc,tauprobe);%Probe pulse with temporal Gaussian profile
xp=x-beta(1);
sigma=(xp*(beta(2))^(-1)).^2;
% Probe=beta(5)*exp(-sigma);%Probe pulse with temporal Gaussian profile
Probe=exp(-sigma);%Probe pulse with temporal Gaussian profile



%Introduction of Heaviside step function in the material response to preserve causality
Th=heaviside(x);%Heaviside step function of t

% Response=beta(5)*Th*(beta(3)-beta(2))^(-1).*(exp(-((x-beta(4))*beta(3)^(-1))-exp(-((x-beta(4))*beta(2)^(-1)))));
% %What is this formula??
% Response=beta(4)*(beta(2)-beta(3))^(-1)*Th.*(-exp(-x*(beta(3)^(-1)))+exp(-x*((beta(2)^(-1)))));
Response=(beta(3)-beta(4))^(-1)*Th.*(-exp(-xp*(beta(4)^(-1)))+exp(-xp*((beta(3)^(-1)))));

% Response=const*Th*(tauth-taur)^(-1).*(exp(-t*tauth^(-1))-exp(-t*taur^(-1)));
% %Ths is the correct formula for the response time


%Taking the Fourier transform
lpad=length(xp);%zero-padding for correct estimate of amplitude in the frequency domain
Ft_Probe=fft(Probe,lpad);%FFT with zero padding
Ft_Probe=Ft_Probe/lpad;%Normalized FFT
Ft_shft_Probe=fftshift(Ft_Probe);%Centered FFT

Ft_Response=fft(Response,lpad);%zero-padding for correct estimate of amplitude in the frequency domain
Ft_Response=Ft_Response/lpad;%Normalized FFT
Ft_shft_Response=fftshift(Ft_Response);%Centered FFT

%Taking the convolution
Mult=Ft_shft_Probe.*Ft_shft_Response;%Taking product of two FT in the frequency domain
Result2=ifft(Mult);%IFFT
Result3=ifftshift(Result2);%Convolution with shifting the center of the outcome. This should be the correct outcome.
Result3=abs(Result3)/max(abs(Result3));
G=1+Result3;

% Result4=conv(Probe,Response,'same');%Direct way of performing convolution in MATLAB without performing FFT
% Result4=Result4/max(Result4);
% G=1+Result4;
end

Answer 1

[Ameer Hamza]

That problem happens because, at some point, the parameters passed to hydrodynamic by fitnlm are such that your function returns NaN or Inf values. When that happens, fitnlm fails and throws that error. You can avoid this by replacing Inf and NaN with very large values, e.g., 1e50. This avoids the issue with Inf while also telling the optimizer that these parameters are not good. Add the following lines at the end of hydrodynamic

G=1+Result3;
G(isinf(G)) = sign(x(isinf(G))).*1e100;
G(isnan(G)) = 1e100;

[Subhajit Bej]

Hi,

I tried that. Though I was able to run the code, the fitting looks weird.

Attached here is the plots. By any chance, this is coming from FFT?

Thanks!

[Ameer Hamza]

It is possibly related to some issue with the hydrodynamic. Somewhere, some values are unexpectedly changing inside the function, which caused a sudden dip in the model. I suggest you add a breakpoint in this function and run it line by line to find out the issue.