Dear all,
I've written a function (attached) to implement a nonparametric
(kernel-based) monotone increasing regression method. Please have a
look and add it to the Statistics package if it seems appropriate.
Best,
Nir
______________
Nir Y. Krakauer
Assistant Professor of Civil Engineering
The City College of New York
New York, NY 10031
(212) 650-8003
nkraka...@ccny.cuny.edu
http://www-ce.ccny.cuny.edu/nir/
## Copyright (C) 2011 Nir Krakauer
##
## This program is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program. If not, see <http://www.gnu.org/licenses/>.
## -*- texinfo -*-
## @deftypefn {Function File} {@var{yy} =} monotone_smooth (@var{x}, @var{y}, @var{h})
## Produce a smooth monotone increasing approximation to a sampled functional dependence y(x) using a kernel method.
##
## @subheading Arguments
##
## @itemize @bullet
## @item
## @var{x} is a vector of values of the independent variable.
##
## @item
## @var{y} is a vector of values of the dependent variable, of the same size as @var{x}. For best performance, it is recommended that the @var{y} already be fairly smooth, e.g. by applying a kernel smoothing to the original values if they are noisy.
##
## @item
## @var{h} is the kernel bandwidth to use. If @var{h} is not given, a "reasonable" value is computed.
##
## @end itemize
##
## @subheading Return values
##
## @itemize @bullet
## @item
## @var{yy} is the vector of smooth monotone increasing function values at @var{x}.
##
## @end itemize
##
## @subheading Examples
##
## @example
## @group
## x = 0:0.1:10;
## y = (x .^ 2) + 3 * randn(size(x)); %typically non-monotonic from the added noise
## ys = ([y(1) y(1:(end-1))] + y + [y(2:end) y(end)])/3; %crudely smoothed via moving average, but still typically non-monotonic
## yy = monotone_smooth(x, ys); %yy is monotone increasing in x
## @end group
##
## @end example
##
## @subheading References
##
## @enumerate
## @item
## Holger Dette, Natalie Neumeyer and Kay F. Pilz (2006), A simple nonparametric estimator of a strictly monotone regression function, @cite{Bernoulli}, 12:469-490
## @item
## Regine Scheder (2007), R Package 'monoProc', Version 1.0-6, @url{http://cran.r-project.org/web/packages/monoProc/monoProc.pdf} (The implementation here is based on the monoProc function mono.1d)
## @end enumerate
## @end deftypefn
## Author: Nir Krakauer <nkraka...@ccny.cuny.edu>
## Description: Nonparametric monotone increasing regression
function yy = monotone_smooth(x, y, h)
n = numel(x);
%set filter bandwidth at a reasonable default value, if not specified
if nargin < 3 || ~isscalar(h)
s = std(x);
h = s / (n^0.2);
end
x_min = min(x);
x_max = max(x);
y_min = min(y);
y_max = max(y);
%transform range of x to [0, 1]
xl = (x - x_min) / (x_max - x_min);
yy = ones(size(y));
%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));
%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
for l = 1:n
tmax = y_max;
tmin = y_min;
wmin = monotone_inverse(tmin);
wmax = monotone_inverse(tmax);
if (wmax == wmin); yy(l) = tmin; end
wt = xl(l);
for i = 1:150
tn = tmin + (wt - wmin) * (tmax -tmin) / (wmax - wmin);
wn = monotone_inverse(tn);
if (abs(wt-wn) < 1E-4) || (tn < (y_min-0.1)) || (tn > (y_max+0.1))
break
end
if wn > wt
tmax = tn;
wmax = wn;
else
tmin = tn;
wmin = wn;
end
end
yy(l) = (tmin + (wt-wmin)*(tmax - tmin)/(wmax- wmin));
end
------------------------------------------------------------------------------
RSA(R) Conference 2012
Save $700 by Nov 18
Register now
http://p.sf.net/sfu/rsa-sfdev2dev1
_______________________________________________
Octave-dev mailing list
Octave-dev@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/octave-dev
