Kernel Bandwidth Optimization 1.0

function [optW, C, W] = sskernel(x,W,str)% [optW, C, W] = sskernel(x,W,str)%% Function `sskernel' returns an 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.%% For more information, visit% http://2000.jukuin.keio.ac.jp/shimazaki/res/kernel.html%% Original paper:% Hideaki Shimazaki and Shigeru Shinomoto% Kernel Bandwidth Optimization in Spike Rate Estimation% Journal of Computational Neuroscience 2010% 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 2010, Hideaki Shimazaki All rights reserved.% http://2000.jukuin.keio.ac.jp/shimazaki%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Parameters Settingsx = 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 FunctionN_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) ));endC = C/2/sqrt(pi);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Optimal Bin Size Selection[optC,nC]=min(C); optW = W(nC);

Requirements:	No special requirements
Platforms:	Matlab
Keyword:	Optimal,  Option,  Optional,  Output,  Reserved,  Resolution,  Sampling,  Smaller,  Specifies,  String
Users rating:	0/10