Custom expression for surface fitting

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!

Answer 1

[Walter Roberson]

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.

Answer 2

[kk]

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?

[Walter Roberson]

The curve fitting routine will only pass in arrays of compatible size. The important point is that the output needs to be exactly the same size as the input.

Can you post your FitFormula function so that we can see how the parameters passed to it by fit() make their way to Streak ?

[kk]

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);
    
while ~isempty(varargin)              %reading optional parameters
    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.