Home > wafo > kdetools > mkernel.m

mkernel

PURPOSE ^

Multivariate Kernel Function.

SYNOPSIS ^

[z,c]=mkernel(varargin)

DESCRIPTION ^

 MKERNEL Multivariate Kernel Function. 
    
  CALL:  z = mkernel(x1,x2,...,xd,kernel); 
         z = mkernel(X,kernel); 
           
  
    z      = kernel function values evaluated at x1,x2,...,xd   
    x1,x2..= input arguments, vectors or matrices with common size 
  or  
    X      = cellarray of vector/matrices with common size 
             (i.e. X{1}=x1, X{2}=x2....) 
  
    kernel = 'epanechnikov'  - Epanechnikov kernel.  
             'epa1'          - product of 1D Epanechnikov kernel.  
             'biweight'      - Bi-weight kernel. 
             'biw1'          - product of 1D Bi-weight kernel. 
             'triweight'     - Tri-weight kernel.   
             'triangular'    - Triangular kernel. 
             'gaussian'      - Gaussian kernel 
             'rectangular'   - Rectangular kernel.  
             'laplace'       - Laplace kernel. 
             'logistic'      - Logistic kernel.   
    
   Note that only the first 4 letters of the kernel name is needed. 
  
  See also  kde, kdefun, kdebin

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function [z,c]=mkernel(varargin) 
002 %MKERNEL Multivariate Kernel Function. 
003 %   
004 % CALL:  z = mkernel(x1,x2,...,xd,kernel); 
005 %        z = mkernel(X,kernel); 
006 %          
007 % 
008 %   z      = kernel function values evaluated at x1,x2,...,xd   
009 %   x1,x2..= input arguments, vectors or matrices with common size 
010 % or  
011 %   X      = cellarray of vector/matrices with common size 
012 %            (i.e. X{1}=x1, X{2}=x2....) 
013 % 
014 %   kernel = 'epanechnikov'  - Epanechnikov kernel.  
015 %            'epa1'          - product of 1D Epanechnikov kernel.  
016 %            'biweight'      - Bi-weight kernel. 
017 %            'biw1'          - product of 1D Bi-weight kernel. 
018 %            'triweight'     - Tri-weight kernel.   
019 %            'triangular'    - Triangular kernel. 
020 %            'gaussian'      - Gaussian kernel 
021 %            'rectangular'   - Rectangular kernel.  
022 %            'laplace'       - Laplace kernel. 
023 %            'logistic'      - Logistic kernel.   
024 %   
025 %  Note that only the first 4 letters of the kernel name is needed. 
026 % 
027 % See also  kde, kdefun, kdebin 
028    
029 %  Reference:   
030 %  B. W. Silverman (1986)  
031 % 'Density estimation for statistics and data analysis'   
032 %  Chapman and Hall, pp. 43, 76  
033 %   
034 %  Wand, M. P. and Jones, M. C. (1995)  
035 % 'Density estimation for statistics and data analysis'   
036 %  Chapman and Hall, pp 31, 103,  175   
037  
038 %Tested on: matlab 5.3 
039 % History: 
040 % Revised pab sep2005 
041 % -replaced reference to kdefft with kdebin 
042 % revised pab aug2005 
043 % -Fixed some bugs 
044 % revised pab Dec2003 
045 % removed some old code   
046 % revised pab 27.04.2001 
047 % - removed some old calls 
048 %  revised pab 01.01.2001 
049 %  - speeded up tri3 
050 %  revised pab 01.12.1999 
051 %   - added four weight, sphere 
052 %   - made comparison smarter => faster execution for d>1 
053 %  revised pab 26.10.1999 
054 %   fixed normalization fault in epan 
055 % by pab 21.09.99   
056 %  added multivariate epan, biweight and triweight 
057 % 
058 % collected all knorm,kepan ... into this file  
059 % adapted from kdetools CB  
060  
061 d=length(varargin)-1; 
062 kstr=varargin{d+1}; % kernel string 
063 if iscell(varargin{1}) 
064   X=varargin{1}; 
065   d=prod(size(X)); 
066 else 
067   X=varargin; 
068 end 
069  
070 switch lower(kstr(1:4))  
071   case {'sphe','epan','biwe','triw','four'} 
072     switch lower(kstr(1:4)) 
073       case 'sphe', r=0;  %Sphere = rect for 1D 
074       case 'epan', r=1;  %Multivariate Epanechnikov kernel. 
075       case 'biwe', r=2;  %Multivariate Bi-weight Kernel 
076       case 'triw', r=3;  %Multi variate Tri-weight Kernel  
077       case 'four', r=4;  %Multi variate Four-weight Kernel 
078         % as r -> infty, b -> infty => kernel -> Gaussian distribution 
079     end  
080     b=1;% radius of the kernel 
081     b2=b^2; 
082     s=X{1}.^2; 
083     k=find(s<=b2); 
084     z=zeros(size(s)); 
085     ix=2; 
086     while (any(k) & (ix<=d)),  
087       s(k)=s(k)+X{ix}(k).^2; 
088       k1=find(s(k)<=b2); 
089       k=k(k1); 
090       ix=ix+1;  
091     end; 
092     if any(k) 
093       c=2^r*prod(1:r)*vsph(d,b)/prod((d+2):2:(d+2*r)); % normalizing constant 
094       %c=beta(r+1,r+1)*vsph(d,b)*(2^(2*r)); % Wand and Jones pp 31 
095       % the commented c above does note yield the right scaling 
096       % for d>1    
097       z(k)=((1-s(k)/b2).^r)/c; 
098     end 
099      
100   case 'rect', % 1D product Rectangular Kernel 
101     z=zeros(size(X{1})); 
102     k=find(abs(X{1})<=1); 
103     ix=2; 
104     while (any(k) & (ix<=d)), 
105       k1 =find(abs(X{ix}(k))<=1); 
106       k=k(k1); 
107       ix=ix+1  
108     end 
109     if any(k)  
110       z(k)=(0.5^d); 
111     end 
112   case {'epa1','biw1','triw1','fou1'} 
113     switch lower(kstr(1:4)) 
114       %case 'rect', r=0;  %rectangular 
115       case 'epa1', r=1;  %1D product Epanechnikov kernel. 
116       case 'biw1', r=2;  %1D product Bi-weight Kernel 
117       case 'tri1', r=3;  %1D product Tri-weight Kernel  
118       case 'fou1', r=4;  %1D product Four-weight Kernel 
119     end  
120     b=1; 
121     b2=b^2; 
122     b21=1/b2; 
123     z=zeros(size(X{1})); 
124     k=find(abs(X{1})<=b); 
125     ix=2;  
126     while (any(k) & (ix<=d)), 
127       %for ix=2:d 
128       k1 =find(abs(X{ix}(k))<=b); 
129       k=k(k1); 
130       ix=ix+1;  
131     end 
132     if any(k) 
133       c=2^r*prod(1:r)*vsph(1,b)/prod((1+2):2:(1+2*r)); % normalizing constant 
134       z(k) = (1-X{1}(k).^2*b21).^r; 
135       for ix=2:d 
136         z(k)=z(k).*(1-X{ix}(k).^2*b21).^r; 
137       end;    
138       z(k)=z(k)/c^d; 
139     end 
140   case 'tria',% 1D product Triangular Kernel  
141     z=zeros(size(X{1})); 
142     k=find(abs(X{1})<1); 
143     ix=2; 
144     while (any(k) & (ix<=d)), 
145       %for ix=2:d 
146       k1 =find(abs(X{ix}(k))<1); 
147       k=k(k1);  
148       ix=ix+1; 
149     end 
150     if any(k)   
151       z(k) = (1-abs(X{1}(k))); 
152       for ix=2:d 
153         z(k)=z(k).*(1-abs(X{ix}(k))); 
154       end 
155     end 
156    case {'norm','gaus'},% multivariate gaussian  Density Function. 
157      s=X{1}.^2; 
158      for ix=2:d   
159        s=s+X{ix}.^2; 
160      end; 
161      z=(2*pi)^(-d/2)*exp(-0.5*s);    
162    case 'lapl' % Laplace Kernel  
163      z=0.5*exp(-abs(X{1})); 
164      for ix=2:d 
165        z=z.*0.5*exp(-abs(X{ix})); 
166      end 
167    case 'logi', % Logistic Kernel  
168      z1=exp(X{1}); 
169      z=z1./(z1+1).^2; 
170      for ix=2:d 
171        z1=exp(X{ix}); 
172        z=z.*z1./(z1+1).^2; 
173      end 
174      
175    otherwise, error('unknown kernel') 
176  end 
177  
178  
179 return

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group


Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003