Kernel Bandwidth Optimization

Detail

Language

English

Title

Kernel Bandwidth Optimization

Free Keywords

bandwidth, density estimation, kernel, mean integrated squared error, optimization, statistics

Description

Function `sskernel' returns optimal bandwidth of the Gauss function in kernel density estimation.

Date

Aug 12, 2009

Last Modified Date

Aug 12, 2009 01:45:55

Created Date

Aug 12, 2009 00:48:31

Contributor

島崎秀昭 (shimazaki)

Item Type

Tool

Change Log(History)

Aug 12, 2009

Modified; Related to.

Tool type

Matlab

Developer

Hideaki Shimazaki

Preview

Tool file

sskernel.m

Type	: text/plain
Size	: 2.7 KB
Last updated	: Aug 12, 2009
Downloads	: 17

Total downloads since Aug 12, 2009 : 17

Readme

function [optW, C, W] = sskernel(x,W,str)
% [optW, C, W] = sskernel(x,W,str)
%
% Function `sskernel' returns optimal bandwidth (standard deviation) 
% of the Gauss density function used in kernel density estimation.
% Optimization principle is to minimize expected L2 loss function between 
% the kernel estimate and an unknown underlying density function. 
% An assumption made is merely that samples are drawn from the density
% independently each other.
%
% The optimal bandwidth is obtained as a minimizer of the formula, 
% sum_{i,j} \int k(x - x_i) k(x - x_j) dx  -  2 sum_{i~=j} k(x_i - x_j), 
% where k(x) is the kernel function. 
%
% Original paper:
% Shimazaki and Shinomoto, Kernel Bandwidth Optimization in Spike Rate Estimation 
% Journal of Computational Neuroscience 2009
% http://dx.doi.org/10.1007/s10827-009-0180-4
%
% Example usage:
% optW = sskernel(x); ksdensity(x,'width',optW); 
%
% Statistics Toolbox is required to execute ksdensity. 
% If it is not available, define the Gauss function as
% `Gauss = @(s,w) 1/sqrt(2*pi)/w*exp(-s.^2/2/w^2);'.
% Computing `mean( Gauss(x-s,optW) )' provides a kernel density estimate at s.
%
% Input argument
% x:    Sample data vector. 
% W (optinal): 
%       A vector of kernel bandwidths.
%       The optimal bandwidth is selected from the elements of W.  
%       Default value is W = logspace(log10(2*dx),log10((x_max - x_min)),50).
%       * Do not search bandwidths smaller than a sampling resolution of data.
% str (optional):	
%       String that specifies the kernel type.
%       This option is reserved for future extention.
%       Default str = 'Gauss'.
%
% Output argument
% optW: Optimal kernel bandwidth.
% W:    Kernel bandwidths examined. 
% C:    Cost functions of W.
%
% See also SSHIST
%
% Copyright (c) 2009, Hideaki Shimazaki All rights reserved.
% http://2000.jukuin.keio.ac.jp/shimazaki

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameters Settings
x = reshape(x,1,numel(x));

str = 'Gauss';

if nargin < 2
    x_min = min(x);
    x_max = max(x);

    buf = sort(abs(diff(sort(x))));
    dx = min(buf(logical(buf ~= 0)));
    Wmin = 2*dx; Wmax = 1*(x_max - x_min);
    W = logspace(log10(Wmin),log10(Wmax),50);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute a Cost Function
N_total = length(x);

tau = triu( ones(N_total,1)*x - x'*ones(1,N_total), 1);
idx = triu( ones(N_total,N_total), 1);
TAU = tau(logical(idx)) .^2;

C = zeros(1,length(W));

for k = 1: length(W)
	w = W(k);
	C(k) = N_total/w + 1/w*sum(sum(2*exp(-TAU/4/w/w) - 4*sqrt(2)*exp(-TAU/2/w/w) )); 
end

C = C/2/sqrt(pi);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Optimal Bin Size Selection
[optC,nC]=min(C); optW = W(nC);

Rights

This work is licensed under a Creative Commons Attribution-NoDerivs 2.5 License

Index

/ Public / Tool

Related to

Item summary

Kernel bandwidth optimization in spike rate estimation.
Shimazaki H , Shinomoto S
Journal of computational neuroscience 2009 [PMID: 19655238]

Histogram Binwidth Optimization
Matlab