changeset 8902:8b328abf65b6 octave-forge

monotone_smooth: check input for extra arguments, return error if h is invalid (instead of silently recalculate it), use printf style for warning message
author carandraug
date Mon, 14 Nov 2011 14:56:10 +0000
parents d93624c61740
children 3c599ff0b35e
files main/statistics/inst/monotone_smooth.m
diffstat 1 files changed, 20 insertions(+), 18 deletions(-) [+]
line wrap: on
line diff
--- a/main/statistics/inst/monotone_smooth.m	Mon Nov 14 14:54:28 2011 +0000
+++ b/main/statistics/inst/monotone_smooth.m	Mon Nov 14 14:56:10 2011 +0000
@@ -75,9 +75,9 @@
 ## Author: Nir Krakauer <nkrakauer@ccny.cuny.edu>
 ## Description: Nonparametric monotone increasing regression
 
-function yy = monotone_smooth(x, y, h)
+function yy = monotone_smooth (x, y, h)
 
-  if nargin < 2 || ~isnumeric(x) || ~isnumeric(y)
+  if nargin < 2 || nargin > 3 || ~isnumeric(x) || ~isnumeric(y)
     print_usage ();
   endif
 
@@ -87,9 +87,11 @@
   x = x(:);
 
   %set filter bandwidth at a reasonable default value, if not specified
-  if nargin < 3 || ~isscalar(h)
+  if (nargin != 3)
     s = std(x);
     h = s / (n^0.2);
+  elseif (!isscalar (h) || !isnumeric (h))
+    error ("third argument 'h' (kernel bandwith) must a numeric scalar")
   end
 
   x_min = min(x);
@@ -106,14 +108,11 @@
   %Epanechnikov smoothing kernel (with finite support)
   %K_epanech_kernel = @(z) (3/4) * ((1 - z).^2) .* (abs(z) < 1);
 
-
-  K_epanech_int = @(z) mean(((abs(z) < 1)/2) - (3/4) * (z .* (abs(z) < 1) - (1/3) * (z.^3) .* (abs(z) < 1)) + (z < -1)); 
+  K_epanech_int = @(z) mean(((abs(z) < 1)/2) - (3/4) * (z .* (abs(z) < 1) - (1/3) * (z.^3) .* (abs(z) < 1)) + (z < -1));
 
   %integral of kernels up to t
   monotone_inverse = @(t) K_epanech_int((y - t) / h);
 
-
-
   %find the value of the monotone smooth function at each point in x
   niter_max = 150; %maximum number of iterations for estimating each value (should not be reached in most cases) 
   for l = 1:n
@@ -129,26 +128,29 @@
       iter_max_reached = 1;
       for i = 1:niter_max
         wt_scaled = (wt - wmin) / (wmax - wmin);
-        tn = tmin + wt_scaled * (tmax - tmin) ;
-        wn = monotone_inverse(tn);
+        tn        = tmin + wt_scaled * (tmax - tmin) ;
+        wn        = monotone_inverse(tn);
         wn_scaled = (wn - wmin) / (wmax - wmin);
-        %if (abs(wt-wn) < 1E-4) || (tn < (y_min-0.1)) || (tn > (y_max+0.1)) %%criterion for break in the R code -- replaced by the following line to hopefully be less dependent on the scale of y
+
+        %if (abs(wt-wn) < 1E-4) || (tn < (y_min-0.1)) || (tn > (y_max+0.1))
+        %% criterion for break in the R code -- replaced by the following line to
+        %% hopefully be less dependent on the scale of y
         if (abs(wt_scaled-wn_scaled) < 1E-4) || (wt_scaled < -0.1) || (wt_scaled > 1.1)
           iter_max_reached = 0;
           break
-        end
+        endif
         if wn > wt
           tmax = tn;
           wmax = wn;
         else
           tmin = tn;
           wmin = wn;
-        end
-      end
+        endif
+      endfor
       if iter_max_reached
-        warning(['at x = ' num2str(x(l)) ', maximum number of iterations reached without convergence; approximation may not be optimal'])
-      end 
+        warning("at x = %g, maximum number of iterations %d reached without convergence; approximation may not be optimal", x(l), niter_max)
+      endif
       yy(l) = tmin + (wt - wmin) * (tmax - tmin) / (wmax - wmin);
-    end
-  end
-end
+    endif
+  endfor
+endfunction