--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/+utils/@math/chi2.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,46 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% DESCRIPTION: Compute chi2 and its gradient.
+%
+% VERSION:     $Id: chi2.m,v 1.1 2011/03/10 16:55:34 congedo Exp $
+%
+% HISTORY:     10-03-2011 G. Congedo
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function [chi2,g] = chi2(p,data,models,Dmodels,lb,ub)
+
+  weights = 1;
+  Nfreeparams = numel(p);
+  Nmdls = numel(models);
+  Ndata = numel(data(:,1));
+  Np = numel(p);
+  
+  mdldata = zeros(Ndata,Nmdls);
+  for kk=1:Nmdls
+    mdldata(:,kk) = models{kk}(p);
+  end  
+  res = (mdldata-data).*weights;
+  if all(p>=lb & p<=ub)
+    chi2 = res'*res;
+    chi2 = sum(diag(chi2));
+  else
+    chi2 = 10e50;
+  end
+  
+  if nargout > 1 % gradient required
+    grad = cell(Nmdls,1);
+    g = zeros(Nmdls,Nfreeparams);
+    for kk=1:Nmdls
+      grad{kk} = zeros(Ndata,Np);
+      for ii=1:Np        
+        grad{kk}(:,ii) = Dmodels{kk}{ii}(p);
+        g(kk,ii) = 2.*res(:,kk)'*grad{kk}(:,ii);
+      end
+    end
+    if Nmdls>1
+      g = sum(g);
+    end
+  end
+  
+end