Trying to create fittype using the convolution of two functions, "Expression is not a valid MATLAB expression, has non-scalar coefficients, or cannot be evaluated"

I am attempting to fit the convolution of a Fermi-Dirac distribution and a Gaussian distribution to some data and am running into issues with defining the fittype inorder to use fit() function. I am able to fit the data using this function when only fitting the Fermi-Dirac distribution but when I attempt to use the convolution in the fittype() function I get an error.

Error using fittype/testCustomModelEvaluation

Expression integral(@(tau)(exp(1.0./x5.^2.*(tau-x7).^2.*(-1.0./2.0)).*(x1+x2./(exp((tau-x3)./(x4.*x6))+1.0)).*3.989422804014327e-1)./x5,-Inf,Inf) is not a valid MATLAB expression, has non-scalar coefficients, or cannot be evaluated:

Bellow is the code I am trying to use along with sample data.

xData = [35.9249999999992	35.9499999999992	35.9749999999992	35.9999999999992	36.0249999999992	36.0499999999992	36.0749999999992	36.0999999999992	36.1249999999992	36.1499999999992	36.1749999999992	36.1999999999992	36.2249999999992	36.2499999999992	36.2749999999992	36.2999999999992	36.3249999999992	36.3499999999992	36.3749999999992	36.3999999999992	36.4249999999992	36.4499999999992	36.4749999999992	36.4999999999992	36.5249999999992	36.5499999999992	36.5749999999992	36.5999999999992	36.6249999999992	36.6499999999992	36.6749999999992	36.6999999999992	36.7249999999992	36.7499999999992	36.7749999999992	36.7999999999992	36.8249999999992	36.8499999999992	36.8749999999992	36.8999999999992	36.9249999999992	36.9499999999992	36.9749999999992	36.9999999999992];
yData = [8547	9125	9066	8723	8455	7687	7797	9793	9825	8438	8030	8007	7433	3543	2352	1573	658	747	945	214	295	185	196	116	109	222	91	373	79	197	137	149	89	142	281	77	159	84	147	287	84	77	95	77];

plot(xData,yData,'.b','HandleVisibility','off')

x = sym('x',[1 7]);

% Function declarations
% Fermi-Dirac
F = x(1)+(x(2)/(exp((x(7)-x(3))/(x(6)*x(4)))+1));
% Gausian 
G = (1/(x(5)*sqrt(2*pi())))*exp(-(1/2)*((x(7))/(x(5)))^2);

% Convolution
syms tau real;

atau = subs(F, x(7), tau); % transform all to the tau axis
btau = subs(G, x(7), x(7) - tau); % give one the t-tau
H = int(atau * btau, tau, -inf, inf);

% Convert to function handle
s = symvar(H);
j = matlabFunction(H);
Coef = cellstr(string([s(1:5)]));
prob = cellstr(string(s(6)));
ind = cellstr(string(s(end)));

ft = fittype(j,...
    'coefficients',Coef,...
    'dependent',{'y'},'independent',ind,...
    'problem',prob);

Error using fittype/testCustomModelEvaluation (line 12)
Expression integral(@(tau)(exp(1.0./x5.^2.*(tau-x7).^2.*(-1.0./2.0)).*(x1+x2./(exp((tau-x3)./(x4.*x6))+1.0)).*3.989422804014327e-1)./x5,-Inf,Inf) is not a valid MATLAB expression, has non-scalar coefficients, or cannot be evaluated:
Error while trying to evaluate FITTYPE function :

Output of the function must be the same size as the input. If FUN is an array-valued integrand, set the 'ArrayValued' option to true.

Error in fittype>iCreateFittype (line 373)
testCustomModelEvaluation( obj );

Error in fittype (line 330)
obj = iCreateFittype( obj, varargin{:} );

Caused by:
Error using fittype/evaluate>I_ERROR_FCN_ (line 179)
Error while trying to evaluate FITTYPE function :

Output of the function must be the same size as the input. If FUN is an array-valued integrand, set the 'ArrayValued' option to true.

Answer 1

[Paras Gupta]

Hi Ethan,

I understand that you get the error provided in the question when you use the 'fittype' function to fit the convolution of a Fermi-Dirac distribution and a Gaussian distribution.

The error occurs because the function handle 'j' contains the 'integral' function with an array-valued integrand but does not have the 'ArrayValued' property set to 'true'. You can refer to the following link for more information on the 'integral' function and 'ArrayValued' property:

https://www.mathworks.com/help/matlab/ref/integral.html

There does not seem to be an option to control the behavior of the 'matlabfunction' function such that the output function handle 'j' has the 'ArrayValued' property set. However, you can try adding the property programatically by modifying the function handle. The code below shows one way to achieve the same:

xData = [35.9249999999992	35.9499999999992	35.9749999999992	35.9999999999992	36.0249999999992	36.0499999999992	36.0749999999992	36.0999999999992	36.1249999999992	36.1499999999992	36.1749999999992	36.1999999999992	36.2249999999992	36.2499999999992	36.2749999999992	36.2999999999992	36.3249999999992	36.3499999999992	36.3749999999992	36.3999999999992	36.4249999999992	36.4499999999992	36.4749999999992	36.4999999999992	36.5249999999992	36.5499999999992	36.5749999999992	36.5999999999992	36.6249999999992	36.6499999999992	36.6749999999992	36.6999999999992	36.7249999999992	36.7499999999992	36.7749999999992	36.7999999999992	36.8249999999992	36.8499999999992	36.8749999999992	36.8999999999992	36.9249999999992	36.9499999999992	36.9749999999992	36.9999999999992];
yData = [8547	9125	9066	8723	8455	7687	7797	9793	9825	8438	8030	8007	7433	3543	2352	1573	658	747	945	214	295	185	196	116	109	222	91	373	79	197	137	149	89	142	281	77	159	84	147	287	84	77	95	77];

% plot(xData,yData,'.b','HandleVisibility','off')

x = sym('x',[1 7]);

% Function declarations
% Fermi-Dirac
F = x(1)+(x(2)/(exp((x(7)-x(3))/(x(6)*x(4)))+1));
% Gausian 
G = (1/(x(5)*sqrt(2*pi())))*exp(-(1/2)*((x(7))/(x(5)))^2);

% Convolution
syms tau real;

atau = subs(F, x(7), tau); % transform all to the tau axis
btau = subs(G, x(7), x(7) - tau); % give one the t-tau
H = int(atau * btau, tau, -inf, inf);

% Convert to function handle
s = symvar(H);
Coef = cellstr(string([s(1:5)]));
prob = cellstr(string(s(6)));
ind = cellstr(string(s(end)));

j = matlabFunction(H)

j = function_handle with value:
    @(x1,x2,x3,x4,x5,x6,x7)integral(@(tau)(exp(1.0./x5.^2.*(tau-x7).^2.*(-1.0./2.0)).*(x1+x2./(exp((tau-x3)./(x4.*x6))+1.0)).*3.989422804014327e-1)./x5,-Inf,Inf)

% Modifying the function handle by converting to string
jStr = functions(j).function; % string equivalent of the function expresion
lastParenPos = find(jStr == ')', 1, 'last'); % Finding the last parenthesis
modifiedStr = insertBefore(jStr, lastParenPos, "', 'ArrayValued', true"); % Inserting property before parenthesis
modifiedJ = str2func(modifiedStr) % converting string back to function handle

modifiedJ = function_handle with value:
    @(x1,x2,x3,x4,x5,x6,x7)integral(@(tau)(exp(1.0./x5.^2.*(tau-x7).^2.*(-1.0./2.0)).*(x1+x2./(exp((tau-x3)./(x4.*x6))+1.0)).*3.989422804014327e-1)./x5,-Inf,Inf','ArrayValued',true)

ft = fittype(modifiedJ,...
    'coefficients',Coef,...
    'dependent',{'y'},'independent',ind,...
    'problem',prob);

Please refer to the following documentations for more information on the functions used above:

• functions - https://www.mathworks.com/help/matlab/ref/functions.html
• insertBefore - https://www.mathworks.com/help/matlab/ref/insertbefore.html
• str2func - https://www.mathworks.com/help/matlab/ref/str2func.html

Hope this helps.