# Custom expression for surface fitting

8 views (last 30 days)
kk on 13 Mar 2019
Commented: kk on 15 Mar 2019
I am sorry, I don't have any formal Matlab training, so I might just be missing a simple point here.
I have a problem figuring out how to write a function for a custom surface-fit equation. In a simple case, I have this equation , where is a log-normal distribution, is a (known) parabola and . Some of the parameters are known, the fit has 6 coefficients to find (3 for the lognormal distribution, ).
If I simply write an expression (FitFormula) for this equation and use
ft = fittype(FitFormula, 'independent', {'x', 'y'}, 'dependent', 'z','problem',{'sigma' ,'y0'});
[fitresult, gof] = fit( [xData, yData], zData, ft, opts, 'problem',{sigma,y0} );
the fit successfully computes the unknown parameters.
But I would need to write a function for this equation rather than a simple in-line expression, because I want to use some more advanced functions involving a convolution. Here is where I cannot figure out how to do it, because the "Streak" function uses both x and y as a variable. So I made sure that it produces a matrix as output:
function y = Streak(x,tau,tau0,w)
C=max([size(tau),size(tau0)]); %number of columns of output
R=max(size(x)); %number of rows in output
if ((~isscalar(x)) && (isrow(x))) x=x'; end;
if C>1 % check if all inputs are vectors
if ((~isscalar(tau)) && (iscolumn(tau))) tau=tau'; end;
if ((~isscalar(tau0)) && (iscolumn(tau0))) tau0=tau0'; end;
if ((~isscalar(w)) && (iscolumn(w))) w=w'; end;
if isscalar(tau) tau=repmat(tau,1,C); end;
if isscalar(tau0) tau0=repmat(tau0,1,C); end;
if isscalar(w) w=repmat(w,1,C); end;
end;
XX0=repmat(x,1,C)-repmat(tau0,R,1);
TT=repmat(tau,R,1);
GG=repmat(w,R,1);
ExpDec=exp(-XX0./TT);
GGauss=exp(GG.^2./(2*TT.^2));
Erf = erf((XX0-GG.^2./TT)./(sqrt(2)*GG))+1;
y = GGauss.*ExpDec.*Erf;
end
When I use this "Streak" function as the FitFormula, Matlab's
function testCustomModelEvaluation( obj )
produces an error
Error using fittype/testCustomModelEvaluation (line 16)
Custom equations must produce an output vector, matrix, or array that is the same size and shape as the input data. This custom equation fails to meet that requirement
even though the size of the output is OK if I simply evaluate it.
Does anyone know what the correct way to define the "Streak" function is? Thaks a lot!

Walter Roberson on 13 Mar 2019
Notice that if your input x is a row vector then you transpose it and work with that to calculate y. That will not be the same orientation as the input. If the input x is a row vector then the output must be a row vector as well.

kk on 13 Mar 2019
Thank you very much for the quick reaction. I'll try to change the function to work with a row x vector. Just a follow-up question: the reason why I changed the output automatically to a matrix is that my data typically have 480 values for x and 50 values for y and they are 3D data (meaning that every combination of x and y has the corresponding z, not that only selected combinations of [x,y] produce a z value). In this case it makes no sense to calculate a formula such as the one above with two row vectors - they would have to have the same size for such a calculation to work. Does your answer mean that the correct way to write the function is to include the possibility of two input row vectors, but the function should produce an error if they have different sizes?
##### 2 CommentsShowHide 1 older comment
kk on 15 Mar 2019
Thank you for your previous answer once again. You were right, I had to make sure that a combination of two input row vectors of the same size produces a row vector of the same size, because this is exactly what Matlab was checking for.
Just to document the solution to the problem: I used an option to optionally keep the original syntax, but basically what the fitting function needed was an analogue of an inline formula:
function y = Streak(x,tau,tau0,w,varargin)
%PL decay including rise time, formula for a convolution of Exp and IRF
% IRF is a Gaussian with standard deviation w: 1/((sqrt(2*pi)*w)*exp(-(x-x0)^2/(2w^2))
% this Gaussian shape should be the normpdf function in Matlab
%input: tau lifetime, tau0 beginning of PL, w pulse width; amplitude and
%background not included
%% option: AutoVectorize, default off, useful for fitting procedures
% if on: any combination of a scalar and vector should work, if a combination of scalars and vectors is used as input, the scalars are replicated into vectors
options = struct( ...
'AutoVectorize', false);
switch upper(varargin{1})
case 'AUTOVECTORIZE'
options.AutoVectorize = varargin{2};
varargin(1:2) = [];
otherwise
error(['Unexpected option: ' varargin{1}])
end
end
C=max([size(tau),size(tau0)]); %number of columns of output
R=max(size(x)); %number of rows in output
if max(C,R)>1 % at least one vector
if ~options.AutoVectorize % use as input, suitable for fitting functions
X=x; TAU0=tau0; TT=tau; GG=w;
else % auto change the output into a corresponding matrix
if (isrow(x)) x=x'; end;
if C>1 % check if all inputs are vectors
if ((~isscalar(tau)) && (iscolumn(tau))) tau=tau'; end;
if ((~isscalar(tau0)) && (iscolumn(tau0))) tau0=tau0'; end;
if ((~isscalar(w)) && (iscolumn(w))) w=w'; end;
if isscalar(tau) tau=repmat(tau,1,C); end;
if isscalar(tau0) tau0=repmat(tau0,1,C); end;
if isscalar(w) w=repmat(w,1,C); end;
end;
X=repmat(x,1,C);
TAU0=repmat(tau0,R,1);
TT=repmat(tau,R,1);
GG=repmat(w,R,1);
end;
try
XX0=X-TAU0;
ExpDec=exp(-XX0./TT);
GGauss=exp(GG.^2./(2*TT.^2));
Erf = erf((XX0-GG.^2./TT)./(sqrt(2)*GG))+1;
y = 1/2*GGauss.*ExpDec.*Erf;
catch ME
switch ME.identifier
case 'MATLAB:dimagree'
error('Matrix dimensions do not agree. Check the input or use the AutoVectorize option.');
otherwise
rethrow(ME);
end;
end;
end
I defined the log-normal distribution in a similar way and I can use this syntax for fitting:
FitFormula='Streak(x,TauY0+TauA*exp(y/TauW),4.5768e-04*y^2-0.6699*y+3.606486e+02,sigma)*LogNormal(y,AA,Ac,Aw) +y0';
ft = fittype(FitFormula, 'independent', {'x', 'y'}, 'dependent', 'z','problem',{'sigma','y0'});
and
[fitresult, gof] = fit( [xData, yData], zData, ft, opts, 'problem', {6.8,0.0078} );
Now I'll see if I can do the same with a convolution.