multivariate kernel density code

[Michael Creel]

Here's a first stab at multivariate kernel density code. I'm planning
on writing a .oct version to allow fitting large data sets without
using up memory (this version uses kronecker products). I'm also
planning on writing leave one out cross validation code. Any comments
on this code are welcome.
Michael

# Copyright (C) 2006 Michael Creel <michael.creel@...>
#
# 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 2 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, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

# kernel_density: multivariate radial symmetric Epanechnikov kernel density estimator
#
# usage: dens = kernel_density(eval_points, data, window_width)
# or kernel_density() for an example
#
# inputs:
# eval_points: PxK matrix of points at which to calculate the density
# data: NxK matrix of data points
# window_width: positive scalar, the smoothing parameter. The fit
# is more smooth as the window width increases.
# outputs:
# dens: Px1 vector giving the fitted density value at each of the
# P evaluation points
#
# Details: multivariate radial symmetric (spherical) Epanechnikov kernel, with
# bandwidth matrix proportional to the Choleski decomposition
# of the inverse of the covariance matrix of the data
# references:
# Wand, M.P. and Jones, M.C. (1995), 'Kernel smoothing', pg. 105.
# http://www.xplore-stat.de/ebooks/scripts/spm/html/spmhtmlframe73.html

function dens = kernel_density(eval_points, data, window_width)

        # generate some data for example plot, if no data provided
        if (nargin == 0)

                # it is interesting to play with this, it has an important effect on the fit
                window_width = 0.5;

                # generate sample of two R.V.s., X and Y
                # X ~ N(0,1), Y ~ Chi squared(2)
                # X, Y are dependently distributed
                d = randn(500,2);
                data = [d(:,1) sumsq(d,2)];

                # evaluation points are on a grid for plotting
                a = (-3:0.1:3)'; # for the N(0,1)
                b = (0:0.1:6)';  # for the Chi squared(2)
                gridsize = rows(a);
                [grid_x, grid_y] = meshgrid(a, b);
                eval_points = [vec(grid_x) vec(grid_y)];
        endif

        nn = rows(eval_points);
        n = rows(data);
        dim = columns(eval_points);

        # blow it up so as to avoid loops (go buy some memory)
        z1  =  kron(ones(n,1), eval_points);
        z2  =  kron(data, ones(nn,1));

        # inverse bandwidth matrix
        P = chol(inv(cov(data)));

    # squared distances between eval_points and data points
        z  =  (z1 - z2)*P/window_width;
        z  =  sumsq(z, 2);

        # Volume of d-dimensional unit sphere
        c = pi ^ (dim/2) / gamma(dim/2 + 1);

        # compute density
        z  =  reshape(z, nn, n);
        dens = (1/2) / c * (dim + 2) * (1 - z);
        dens = dens .* (z < 1);
    dens  = sum(dens,2) / (n * window_width ^ (dim));
    dens = dens*det(P);

        if (nargin == 0)
                dens = reshape(dens, gridsize, gridsize);
                title("Kernel density fit: dependent N(0,1) and Chi^2(2) data");
                xlabel("true marginal density is N(0,1)");
                ylabel("true marginal density is Chi^2(2)");
                legend("off");
                surf(grid_x, grid_y, dens);
        endif
endfunction

_______________________________________________
Octave-sources mailing list
Octave-sources@...
https://www.cae.wisc.edu/mailman/listinfo/octave-sources