density.m (915B)
1 function f=density(X,t,w,h); 2 3 % Kernel estimation of the density, based on conversation w/ Frank Wolak. 4 % 5 % X == Data, Nx1 6 % t == range over which to compute the density, e.g. -1 to +1, step .1 7 % w == weights (optional) s.t. sum(w)=1 e.g. w=1/N 8 % h == bandwidth. If h=0 or doesn't exist, then use 9 % 10 % h=0.9*std(X)*N^(-1/5) 11 % 12 % The density is constructed as 13 % 14 % f(t) = 1/h* SUM(i=1:N) w_i K( (t-X_i)/h ) 15 % 16 % where K(.) is the kernel, assumed to be a standard normal. 17 18 N=length(X); 19 if exist('h')~=1; 20 h=0.9*std(X)*N^(-1/5); 21 end; 22 if exist('w')~=1; 23 w=ones(N,1)/N; 24 end; 25 fprintf('Bandwidth h = %10.5f\n',h); 26 27 T=length(t); 28 f=zeros(T,1); 29 30 for j=1:T; 31 xpoint=(t(j)-X)/h; 32 f(j)=1/h*sum(w.*normpdf(xpoint)); 33 end; 34 35 % There's a slight issue here about what the area equals... 36 fprintf('Sum of f without normalizing: %8.4f\n',sum(f)); 37 %%f=f/sum(f); % Normalize so that it adds to one