Learn how easy it is to sync an existing GitHub or Google Code repo to a SourceForge project! See Demo

Close

Diff of /inst/rgdtsmcore.m [96593a] .. [ed273d] Maximize Restore

  Switch to side-by-side view

--- a/inst/rgdtsmcore.m
+++ b/inst/rgdtsmcore.m
@@ -36,8 +36,11 @@
 ##  two element vector [min,max] giving the desired output range of
 ##  @var{xhat}; if not provided or if the provided range is less than
 ##  the range of the data, the default is the min and max of @var{x}
+## @item "weights", @var{vector}
+##  A vector of weighting values for fitting each point in the data.
 ## @item "relative"
-##  use relative differences for the goodnes of fit term
+##  use relative differences for the goodnes of fit term.  Conflicts
+##  with the "weights" option.
 ## @item "midpointrule"
 ##  use the midpoint rule for the integration terms rather than a direct sum
 ##@end table
@@ -54,16 +57,10 @@
   Nhat = 100;
   range = [min(x),max(x)];
   gridx = 0;
+  weights = 0;
   relative = 0;
   midpr = 0;
   
-  ## this is a hack to allow passing varargin from another function...
-  ## is there a better way?
-  if length(varargin)
-    if iscell(varargin{1})
-      varargin = varargin{1};
-    endif
-  endif
   ## parse the provided options
   if ( length(varargin) )
     for i = 1:length(varargin)
@@ -78,6 +75,9 @@
 	  case "range"
 	    gridx = 1;
 	    range = varargin{i+1};
+	  case "weights"
+	    weights = 1;
+	    weightv = varargin{i+1};
 	  case "relative"
 	    relative = 1;
 	  case "midpointrule"
@@ -91,6 +91,9 @@
   if (gridx & midpr)
     warning("midpointrule is not used if mapping to regular spaced x")
     midpr = 0;
+  endif
+  if (weights & relative)
+    warning("relative differences is used if a weighting vector is provided")
   endif
   
   N = length(x);
@@ -115,8 +118,11 @@
   endif
 
   ## construct "weighting" matrices W and U
-  ## use relative differences
-  if (relative)
+  if (weights)
+    ## use arbitrary weighting as provided
+    W = diag(weightv);
+  elseif (relative)
+    ## use relative differences
     Yinv = spdiag(1./y);
     W = Yinv^2;
   else
@@ -143,15 +149,18 @@
   endif
   
   ## Smoothing
-  delta = trace(D'*D)/Nhat^(2+d);  # using "relative" affects this!
-  ##yhat = (M'*W*M + lambda*D'*U*D) \ M'*W*y;
-  C = chol(M'*W*M + lambda*delta^(-1)*D'*U*D);
-  yhat = C \ (C' \ M'*W*y);
+  delta = trace(D'*D)/Nhat^(2+d);  # using "relative" or other weighting affects this!
+  yhat = (M'*W*M + lambda*delta^(-1)*D'*U*D) \ M'*W*y;
+  #[R,P,S] = splchol(M'*W*M + lambda*delta^(-1)*D'*U*D);
+  #yhat = S*(R'\(R\(S'*M'*W*y)));
   
   ## Computation of hat diagonal and cross-validation
   if (nargout > 2)
     ## from AIChE J. (2006) 52, 325
-    H = M * ( C \ (C' \ M'*W) );
+    ## note: chol factorization does not help speed up the computation of H;
+    ## should implement Eiler's partial H computation if many point smoothing by GCV is needed
+    ##H = M*(S*(R'\(R\(S'*M'*W))));
+    H = M*((M'*W*M + lambda*delta^(-1)*D'*U*D)\M'*W);
     ## note: this is variance, squared of the standard error that Eilers uses
     v = (M*yhat - y)'*(M*yhat - y)/N / (1 - trace(H)/N)^2;
   endif