Adding multidimensional gaussian fitting tool from FEX

fitting/gaussfitn.m

@@ -0,0 +1,224 @@
+function [params,varargout] = gaussfitn(xdata,zdata,params0, lbCell, ubCell,varargin)
+% N-dimensional fitting of a Gaussian+constant model from scattered data.
+%
+%This function uses lsqcurvefit to fit parameters D, A, mu, sig to the R^N-->R
+%model function,
+%
+%        z(x) = D + A*exp( -0.5 * (x-mu).' * inv(sig) *(x-mu) )
+%
+%Here A and D are unknown scalars, mu is an unknown Nx1 mean vector, and sig is an
+%unknown NxN covariance matrix.
+%
+%SYNTAX:
+%
+%  [params,resnorm, residual,exitflag,output] = gaussfitn(xdata,zdata,params0,LB,UB,Name,Value)
+%              
+%INPUTS (required):           
+%              
+%      xdata:  MxN matrix whose rows specify M scattered samples in R^N
+%
+%      zdata:  Mx1 vector of corresponding samples z(xdata) 
+%
+%INPUTS (optional)
+%
+%    params0:  Cell array of initial parameter estimates {D0,A0,mu0,sig0}.
+%              Can also be empty [] in which case default initial guesses
+%              are autogenerated. Can also consist of cell array of empty
+%              and non-empty elements like {D0,[],mu0,[]} in which case
+%              default initial guesses are generated for select parameters.
+%              
+%         LB:  Cell array of lower bounds {D_LB, A_LB, mu_LB} on D, A, and mu.
+%
+%         UB:  Cell array of upper bounds {D_UB, A_UB, mu_UB} on D, A, and mu.       
+%
+% Name,Value:  Name/Value option pairs compatible with lsqcurvefit. See,
+%              <https://www.mathworks.com/help/optim/ug/lsqcurvefit.html#buuhcjo-options>.
+%              By default, however, SpecifyConstraintGradient=true unless
+%              over-ridden.
+%              
+%              
+%OUTPUTS:          
+%              
+%           A: Final estimate of the parameters as a cell array {D,A,mu,sig}
+%     resnorm: As in lsqcurvefit
+%    residual: As in lsqcurvefit
+%    exitflag: As in lsqcurvefit
+%      output: As in lsqcurvefit
+   if nargin<3 || isempty(params0)
+       params0={[],[],[],[]};
+   end
+
+   if nargin<4
+      lbCell=[]; 
+   end
+
+   if nargin<5
+      ubCell=[]; 
+   end
+   [M,N]=size(xdata);
+
+   validateattributes(xdata,{'numeric'},{'nonempty'},'','xdata');
+   validateattributes(zdata,{'numeric'},{'numel',M},'','zdata');
+
+   zdata=zdata(:);
+
+    %% Parse initial parameter guesses
+
+    [D0,A0,mu0,sig0]=deal(params0{:});
+
+    wts=zdata-min(zdata); wts=wts/sum(wts);
+
+    if isempty(mu0)
+       mu0=wts.'*xdata; 
+    end
+    mu0=mu0(:);
+    validateattributes(mu0,{'numeric'},  {'numel',N},'','mu0');
+
+
+    if isempty(sig0)
+        dx=xdata-mu0.';
+        sig0=inv(dx.'*(wts.*dx));
+       C0=chol(sig0,'lower');
+    else
+       validateattributes(sig0,{'numeric'},  {'square','nrows',N},'','sig0');
+       sig0=(sig0+sig0.')/2;
+       C0= chol(inv(sig0),'lower');
+    end
+    Lmap=tril(true(N));
+    L0=C0(Lmap);
+
+
+    p0=[0;1;mu0;L0];
+
+    if isempty([D0,A0]) %both empty
+        ExpTerms=objfun(p0,xdata,Lmap);
+        lsol=(ExpTerms.^[0,1])\zdata;
+        D0=lsol(1); A0=lsol(2);
+    elseif isempty(A0) %but D0 non-empty
+        ExpTerms=objfun(p0,xdata,Lmap);
+        A0=ExpTerms\(zdata-D0);
+    elseif isempty(D0) %but A0 non-empty
+        ExpTerms=objfun(p0,xdata,Lmap);
+        D0=mean(A0*ExpTerms);
+    end
+
+    validateattributes(D0,{'numeric'},{'scalar'},'','D0');
+    validateattributes(A0,{'numeric'},{'scalar'},'','A0');
+
+    p0(1:2)=[D0,A0];
+
+    %% Parse upper and lower bounds
+
+
+    if iscell(lbCell) && numel(lbCell)<4
+       lbCell{4}=[];
+    elseif isempty(lbCell)
+       lbCell={[],[],[],[]}; 
+    end
+
+    if iscell(ubCell) && numel(ubCell)<4
+       ubCell{4}=[];    
+    elseif isempty(ubCell)
+       ubCell={[],[],[],[]}; 
+    end   
+
+    lb=dealBounds('LB',Lmap, lbCell{:}); %support bounds on L, though undocumented
+    ub=dealBounds('UB',Lmap, ubCell{:});
+
+
+
+    %% Do the optimization
+
+    opts=optimoptions(@lsqcurvefit,'SpecifyObjectiveGradient',true);
+    opts=optimoptions(opts,varargin{:});
+
+   [p,varargout{1:nargout-1}] = ... %lambda nad Jacobian can also be returned (undocumented)
+       lsqcurvefit(@(p,xd) objfun(p,xd,Lmap),p0,xdata,double(zdata),lb,ub,opts);  
+    D=p(1);
+    A=p(2);
+    mu=p(3:N+2);
+    L=zeros(N);
+     L(Lmap)=p(N+3:end);
+
+    sig=inv(L*L.');
+
+    params={D,A,mu,sig};
+
+
+    function [z,Jac] = objfun(p,xdata,Lmap)
+   [M,N]=size(xdata);
+    p=p(:);
+
+    D=p(1);
+    A=p(2);
+
+    mu=p(3:N+2);
+
+    C=zeros(N);
+     C(Lmap)=p(N+3:end); 
+
+    Q=C*C.'; 
+
+   dx = (xdata-mu.');
+   dxQ = dx*Q;
+
+
+   Jac=cell(1,4);
+   Jac{2} = exp( -sum( dxQ.*dx , 2)/2 );
+
+   ze=A*Jac{2};
+   z=ze+D;
+
+
+   if nargout<2, return; end
+
+   %% Jacobian computation
+
+   Jac{1}=ones(M,1);
+
+   Jac{3}=dxQ.*ze;
+
+   tmp=reshape( dx*C,[],1,N).*dx;
+   tmp=reshape(tmp,[],N^2);
+
+   Jac{4}= (-ze).*tmp(:,Lmap(:));
+
+   Jac=cell2mat(Jac);
+
+
+
+ function bound_vector=dealBounds(bType,Lmap, D, A, mu, C)  
+
+     switch bType
+
+         case 'LB'
+
+             sign=-1;
+
+         case 'UB'
+
+             sign=+1;
+
+         otherwise
+
+             error 'Unrecognized'
+
+     end
+
+
+     N=size(Lmap,1);
+
+     if isempty(D), D=sign*inf; end
+     if isempty(A), A=sign*inf; end    
+     if isempty(mu), mu=sign*inf(N,1); end
+     if isempty(C), C=sign*inf(N); end        
+
+     validateattributes(D,{'numeric'}, {'scalar'},'',['D_' bType]);
+     validateattributes(A,{'numeric'}, {'scalar'},'',['A_' bType]);
+     validateattributes(mu,{'numeric'},{'numel',N},'',['mu_' bType]);
+     validateattributes(C,{'numeric'}, {'numel',N^2},'',['C_' bType]);
+
+     L=C(Lmap);
+
+     bound_vector=[D;A;mu(:);L(:)];
+