X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?p=CreaPhase.git;a=blobdiff_plain;f=octave_packages%2Feconometrics-1.0.8%2Fkernel_example.m;fp=octave_packages%2Feconometrics-1.0.8%2Fkernel_example.m;h=20bf39947a7b63bfdd4e588c8ad0aeac3e78c6c6;hp=0000000000000000000000000000000000000000;hb=f5f7a74bd8a4900f0b797da6783be80e11a68d86;hpb=1705066eceaaea976f010f669ce8e972f3734b05 diff --git a/octave_packages/econometrics-1.0.8/kernel_example.m b/octave_packages/econometrics-1.0.8/kernel_example.m new file mode 100644 index 0000000..20bf399 --- /dev/null +++ b/octave_packages/econometrics-1.0.8/kernel_example.m @@ -0,0 +1,140 @@ +# Copyright (C) 2007 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, see . + +# kernel_example: examples of how to use kernel density and regression functions +# requires the optim and plot packages from Octave Forge +# +# usage: kernel_example; + +# sample size (default n = 500) - should get better fit (on average) +# as this increases, supposing that you allow the optimal window width +# to be found, by uncommenting the relevant lines +n = 500; + +# set this to greater than 0 to try parallel computations (requires MPITB) +compute_nodes = 0; +nodes = compute_nodes + 1; # count master node + +close all; +hold off; + +############################################################ +# kernel regression example +# uniformly spaced data points on [0,2] +x = 1:n; +x = x'; +x = 2*x/n; +# generate dependent variable +trueline = x + (x.^2)/2 - 3.1*(x.^3)/3 + 1.2*(x.^4)/4; +sig = 0.5; +y = trueline + sig*randn(n,1); +tic; +fit = kernel_regression(x, y, x); # use the default bandwidth +t1 = toc; +printf("\n"); +printf("########################################################################\n"); +printf("time for kernel regression example using %d data points and %d compute nodes: %f\n", n, nodes, t1); +plot(x, fit, ";fit;", x, trueline,";true;"); +grid("on"); +title("Example 1: Kernel regression fit"); + +############################################################ +# kernel density example: univariate - fit to Chi^2(3) data +data = sumsq(randn(n,3),2); +# evaluation point are on a grid for plotting +stepsize = 0.2; +grid_x = (-1:stepsize:11)'; +bandwidth = 0.55; +# get optimal bandwidth (time consuming, uncomment if you want to try it) +# bandwidth = kernel_optimal_bandwidth(data); +# get the fitted density and do a plot +tic; +dens = kernel_density(grid_x, data, bandwidth, "__kernel_normal", false, false, compute_nodes); +t1 = toc; +printf("\n"); +printf("########################################################################\n"); +printf("time for univariate kernel density example using %d data points and %d compute nodes: %f\n", n, nodes, t1); +printf("A rough integration under the fitted univariate density is %f\n", sum(dens)*stepsize); +figure(); +plot(grid_x, dens, ";fitted density;", grid_x, chi2pdf(grid_x,3), ";true density;"); +title("Example 2: Kernel density fit: Univariate Chi^2(3) data"); + +############################################################ +# kernel density example: bivariate +# X ~ N(0,1) +# Y ~ Chi squared(3) +# X, Y are dependent +d = randn(n,3); +data = [d(:,1) sumsq(d,2)]; +# evaluation points are on a grid for plotting +stepsize = 0.2; +a = (-5:stepsize:5)'; # for the N(0,1) +b = (-1:stepsize:9)'; # for the Chi squared(3) +gridsize = rows(a); +[grid_x, grid_y] = meshgrid(a, b); +eval_points = [vec(grid_x) vec(grid_y)]; +bandwidth = 0.85; +# get optimal bandwidth (time consuming, uncomment if you want to try it) +# bandwidth = kernel_optimal_bandwidth(data); +# get the fitted density and do a plot +tic; +dens = kernel_density(eval_points, data, bandwidth, "__kernel_normal", false, false, compute_nodes); +t1 = toc; +printf("\n"); +printf("########################################################################\n"); +printf("time for multivariate kernel density example using %d data points and %d compute nodes: %f\n", n, nodes, t1); +dens = reshape(dens, gridsize, gridsize); +printf("A rough integration under the fitted bivariate density is %f\n", sum(sum(dens))*stepsize^2); +figure(); +surf(grid_x, grid_y, dens); +title("Example 3: Kernel density fit: dependent bivariate data"); +xlabel("true marginal density is N(0,1)"); +ylabel("true marginal density is Chi^2(3)"); +# more extensive test of parallel +if compute_nodes > 0 # only try this if parallel is available + ns =[4000; 8000; 10000; 12000; 16000; 20000]; + printf("\n"); + printf("########################################################################\n"); + printf("kernel regression example with several sample sizes serial/parallel timings\n"); + figure(); + clg; + title("Compute time versus nodes, kernel regression with different sample sizes"); + xlabel("nodes"); + ylabel("time (sec)"); + hold on; + ts = zeros(rows(ns),4); + for i = 1:rows(ns) + for compute_nodes = 0:3 + nodes = compute_nodes + 1; + n = ns(i,:); + x = 1:n; + x = x'; + x = 2*x/n; + # generate dependent variable + trueline = x + (x.^2)/2 - 3.1*(x.^3)/3 + 1.2*(x.^4)/4; + sig = 0.5; + y = trueline + sig*randn(n,1); + bandwidth = 0.45; + tic; + fit = kernel_regression(x, y, x, bandwidth, "__kernel_normal", false, false, compute_nodes); + t1 = toc; + ts(i, nodes) = t1; + plot(nodes, t1, "*"); + printf(" %d data points and %d compute nodes: %f\n", n, nodes, t1); + endfor + plot(ts(i,:)'); + endfor + hold off; +endif