how to fit implicit function to data

I have these data points

data=[0   1
    0.355257 0.726061
    0.749002 0.382696
    1        0
    2.401699 -0.9817
    3.049682 -1.49588
    4.105308 -2.51708
    3.419108 -2.79514
    1.784289 -2.188
    0.763907 -1.56124
    0        -1];

I need to fit this implicit function to the data

f = @(x,y) abs(x)^a + abs(y)^b +c*x*y-1;

where a, b, c are the 3 positive parameters to be estimated.

Can anyone help me?

[Torsten]

Use lsqnonlin.

Answer 1

[John D'Errico]

I'll just throw it into the curve fitting toolbox, although any tool would have sufficed for a first try. You should recognize that other tools, such as lsqnonlin from the optimization toolbox, nlinfit from stats TB, etc., would not require a major change to make them solve this too.

z = ones(11,1);
ft = fittype( @(a,b,c,x,y) abs(x).^a + abs(y).^b +c*x.*y, ...
            'independent', {'x', 'y'},'dependent', 'z' );
mdl = fit(data,z,ft,'startpoint',[1 1 1])
     General model:
     mdl(x,y) = abs(x).^a+abs(y).^b+c*x.*y
     Coefficients (with 95% confidence bounds):
       a =       1.868  (1.805, 1.932)
       b =       2.355  (2.27, 2.441)
       c =       2.112  (1.944, 2.279)
       
mdl(data)
ans =
            1
       1.1598
       1.2921
            1
       1.1184
      0.97943
      0.96971
       1.0215
       1.0299
      0.94193
            1

It seems to have worked reasonably. Note the CFT does not provide for bound constraints on the parameters. So if any of a,b,c were going negative, I would have needed to use another tool for the task, or been more creative in my use of fit.

One nice feature of the CFT is it tells me if the parameters are going near zero, in the case of parameters that are desired to bo NOT less than zero. Here, the confidence bounds provided by fit tell me that it would not be considering zero for any of them.

A problem with a fit like this, is it is not really what you probably want, in the sense that it would be nice if it were some sort of orthogonal fit, minimizing the an orthogonal distance of the points to the curve. As well, the fitting tool makes implicit assumptions about the "noise" structrue in the data, that are surely not valid in this case. Remember that z (the target here) was constant at 1. What noise? Don't got no steenkin' noise. Not in z at least.

So, really this is a variation of an errors in variables problem. And one thing you should know about errors in variables problems, is they can be nasty as hell to solve. As well, the result can have biases in the parameters, that there will be specific points that will have considerably more weight in your final result than others, and noise in your data (thus in x and in y) can do strange things.

Regardless, did it work?

syms X Y
fimplicit(abs(X).^mdl.a + abs(Y).^mdl.b + mdl.c*X.*Y - 1)
hold on
grid on
plot(x,y,'ro')

Could be better, could be worse. As I said though, true orthogonal fits are way more difficult to implement than the simple thing I did here.

[Alessandraro]

thank you very much

Answer 2

[Rik]

Or, if you don't have the curve fitting toolbox (or the optimization toolbox), you can use fminsearch. The code below might be a bit overpowered for your use case, but it is re-usable. Both functions shown below are meant to end up on the FEX at some point.

data=[0   1
    0.355257 0.726061
    0.749002 0.382696
    1        0
    2.401699 -0.9817
    3.049682 -1.49588
    4.105308 -2.51708
    3.419108 -2.79514
    1.784289 -2.188
    0.763907 -1.56124
    0        -1];
x=data(:,1);y=data(:,2);
f = @(fitval,x) abs(x).^fitval(1) + abs(y).^fitval(2) +fitval(3).*x.*y-1;
initial_guess=[1 1 1];
bounds=[0 inf;0 inf;0 inf];
fit_val=fit_fminsearch(f,x,y,initial_guess,bounds);

clc
fprintf('a=%.4f\nb=%.4f\nc=%.4f\n',fit_val)

function [fit_val,gof_struct]=fit_fminsearch(f,x,y,initial_guess,bounds)
%Use fminsearch to fit data to a function
%The first input should be a function handle, anonymous function or inline
%function (for releases without anonymous functions) which accepts the fit
%value as the first input and the x as the second input. It must allow x
%and initial_guess as inputs and return a vector the same size as y.
%
%The bounds input is optional and should be a matrix with two columns and
%length(initial_guess) rows. By default this is set to -inf to inf.
%
% Examples:
%
% f=inline('b*x','b','x');
% x=1:10;y=2*x;initial_guess=1;
% B=fit_fminsearch(f,x,y,initial_guess);
% 
% f=@(b,x) b*x;
% x=1:10;y=2*x;initial_guess=1;
% B=fit_fminsearch(f,x,y,initial_guess);

%#dependencies{ifversion}
persistent OLS
if isempty(OLS)
    % adapted Ordinary Least Squares cost function
    OLS_as_char=[...
        'sum((f(b,x) - y).^2) +'... % OLS part
        'inf^(any(b<bounds(:,1))-1) +'...%converts true to inf
        'inf^(any(b>bounds(:,2))-1)'];
    if ifversion('<',7,'Octave','<','1')
        %Anonymous functions were not introduced, so use an inline function
        %instead. This will be removed in later releases, so avoid inline
        %when this is possible.
        OLS=inline(OLS_as_char,'b','x','y','f','bounds'); %#ok<DINLN>
    else
        %The syntax checker of older releases trip up over statements
        %defining anonymous function, so use eval to trick it.
        eval(['OLS=@(b,x,y,f,bounds) ' OLS_as_char ';'])
    end
end

if nargin==4
    %set the default bounds
    bounds=inf*ones(numel(initial_guess),2);
    bounds(:,1)=-inf;
end

opts = optimset('MaxFunEvals',50000, 'MaxIter',10000);
% Use 'fminsearch' to minimise the 'OLS' function
fit_val = fminsearch(OLS, initial_guess(:), opts,x,y,f,bounds);

if nargout==2
    err=y-f(fit_val,x);
    SSres=sum(err.^2);
    SStot=sum((y-mean(y)).^2);
    rsq=1-(SSres./SStot);
    RMS=sqrt(mean(err.^2));
    
    gof_struct=struct;
    gof_struct.sse=SSres;
    gof_struct.rsquare=rsq;
    gof_struct.rmse=RMS;
end
end
function tf=ifversion(test,Rxxxxab,varargin)
%Determine if the current version satifies a version restriction
%
% To keep the function fast, no input checking is done.
%
% Syntax:
% tf=ifversion(test,Rxxxxab)
% tf=ifversion(test,Rxxxxab,'Octave',test_for_Octave,v_Octave)
%
% Output:
% tf       - If the current version satisfies the test this returns true.
%            This works similar to verLessThan.
%
% Inputs:
% Rxxxxab - Char array containing a release description (e.g. 'R13',
%           'R14SP2' or 'R2018b') or the numeric version.
% test    - Char array containing a logical test. The interpretation of
%           this is equivalent to eval([current test Rxxxxab]). For
%           examples, see below.
%
% Examples:
% ifversion('>=','R2009a') returns true when run on R2009a or later
% ifversion('<','R2016a') returns true when run on R2015b or older
% ifversion('==','R2018a') returns true only when run on R2018a
% ifversion('<',0,'Octave','>',0) returns true only on Octave
%
% The conversion is based on a manual list and therefore needs to be
% updated manually, so it might not be complete. Although it should be
% possible to load the list from Wikipedia, this is not implemented.
%
% Version: 0.1
% Date:    2018-10-27
% Author:  H.J. Wisselink
% Licence: CC by-nc-sa 4.0 ( creativecommons.org/licenses/by-nc-sa/4.0 )
% Email=  'h_j_wisselink*alumnus_utwente_nl';
% Real_email = regexprep(Email,{'*','_'},{'@','.'})

%The decimal of the version numbers are padded with a 0 to make sure v7.10
%is larger than v7.9. This does mean that any numeric version input needs
%to be adapted.
persistent v_num v_dict octave
if isempty(v_num)
    %test if Octave is used instead of Matlab
    octave=exist('OCTAVE_VERSION', 'builtin');
    
    %get current version number
    v_num=version;
    ind=strfind(v_num,'.');
    if numel(ind)~=1,v_num(ind(2):end)='';end
    v_num=[str2double(v_num(1:(ind-1))) str2double(v_num((ind+1):end))];
    v_num=v_num(1)+v_num(2)/100;
    
    %get dictionary to use for ismember
    v_dict={'R13' 6.05;'R13SP1' 6.05;'R13SP2' 6.05;'R14' 7;'R14SP1' 7;...
        'R14SP2' 7.00;'R14SP3' 7.01;'R2006a' 7.02;'R2006b' 7.03;...
        'R2007a' 7.04;'R2007b' 7.05;'R2008a' 7.06;'R2008b' 7.07;...
        'R2009a' 7.08;'R2009b' 7.09;'R2010a' 7.10;'R2010b' 7.11;...
        'R2011a' 7.12;'R2011b' 7.13;'R2012a' 7.14;'R2012b' 8.00;...
        'R2013a' 8.01;'R2013b' 8.02;'R2014a' 8.03;'R2014b' 8.04;...
        'R2015a' 8.05;'R2015b' 8.06;'R2016a' 9.00;'R2016b' 9.01;...
        'R2017a' 9.02;'R2017b' 9.03;'R2018a' 9.04;'R2018b' 9.05};
end

if octave
    if nargin==2
        warning('HJW:ifversion:NoOctaveTest',['No version test ',...
            'for Octave was provided.' char(10) 'This function ',...
            'might return an unexpected outcome.']) %#ok<CHARTEN>
        %Use the same test as for Matlab, which will probably fail.
        L=ismember(v_dict(:,1),Rxxxxab);
        if sum(L)~=1
            v=NaN;
            warning('HJW:ifversion:NotInDict',...
                'The requested version is not in the hard-coded list.')
        else
            v=v_dict{L,2};
        end
    else
        [test,v]=deal(varargin{4:5});
        %convert 4.1 to 4.01
        v=0.1*v+0.9*fix(v);
    end
else
    if isnumeric(Rxxxxab)
        %convert R notation to numeric and convert 9.1 to 9.01
        v=0.1*Rxxxxab+0.9*fix(Rxxxxab);
    else
        L=ismember(v_dict(:,1),Rxxxxab);
        if sum(L)~=1
            v=NaN;
            warning('HJW:ifversion:NotInDict',...
                'The requested version is not in the hard-coded list.')
        else
            v=v_dict{L,2};
        end
    end
end
switch test
    case '=='
        tf= v_num == v;
    case '<'
        tf= v_num <  v;
    case '<='
        tf= v_num <= v;
    case '>'
        tf= v_num >  v;
    case '>='
        tf= v_num >= v;
end
end