From 4e17882dd5fadc94e7b668f33a27aedb9f04e882 Mon Sep 17 00:00:00 2001 From: Jan Benda Date: Mon, 27 Nov 2017 10:26:03 +0100 Subject: [PATCH] some code fixes --- statistics/code/cumulative.out | 6 +++++ statistics/code/gaussiankerneldensity.m | 36 ++++++++++++++----------- statistics/lecture/kerneldensity.py | 2 +- 3 files changed, 28 insertions(+), 16 deletions(-) create mode 100644 statistics/code/cumulative.out diff --git a/statistics/code/cumulative.out b/statistics/code/cumulative.out new file mode 100644 index 0000000..ed8ab1d --- /dev/null +++ b/statistics/code/cumulative.out @@ -0,0 +1,6 @@ +>> cumulative +data : probability of x<-1: 0.15 +gauss: probability of x<-1: 0.16 + +data : 5% percentile at -1.77 +gauss: 5% percentile at -1.65 \ No newline at end of file diff --git a/statistics/code/gaussiankerneldensity.m b/statistics/code/gaussiankerneldensity.m index 4cb7925..8f97818 100644 --- a/statistics/code/gaussiankerneldensity.m +++ b/statistics/code/gaussiankerneldensity.m @@ -1,15 +1,15 @@ -data = randn(100, 1); % generate some data -sigma = 0.2; % standard deviation of Gaussian kernel -xmin = -4.0; % minimum x value for kernel density -xmax = 4.0; % maximum x value for kernel density -dx = 0.05*sigma; % step size for kernel density -xg = [-4.0*sigma:dx:4.0*sigma]; % x-axis for single Gaussian kernel +data = randn(100, 1); % generate some data +sigma = 0.2; % std. dev. of Gaussian kernel +xmin = -4.0; % minimum x value for kernel density +xmax = 4.0; % maximum x value for kernel density +dx = 0.05*sigma; % step size for kernel density +xg = [-4.0*sigma:dx:4.0*sigma]; % x-axis for single Gaussian kernel % single Gaussian kernel: kernel = exp(-0.5*(xg/sigma).^2)/sqrt(2.0*pi)/sigma; -ng = (length(kernel)-1)/2; % half the length of the Gaussian -x = [xmin:dx:xmax+0.5*dx]; % x-axis for kernel density -kd = zeros(1, length(x)); % vector for kernel density -for i = 1:length(data) % for every data value ... +ng = floor((length(kernel)-1)/2); % half the length of the Gaussian +x = [xmin:dx:xmax+0.5*dx]; % x-axis for kernel density +kd = zeros(1, length(x)); % vector for kernel density +for i = 1:length(data) % for every data value ... xd = data(i); % index of data value in kernel density vector: inx = round((xd-xmin)/dx)+1; @@ -17,13 +17,13 @@ for i = 1:length(data) % for every data value ... k0 = inx-ng; % end index for Gaussian in kernel density vector: k1 = inx+ng; - g0 = 1; % start index in Gaussian + g0 = 1; % start index in Gaussian g1 = length(kernel); % end index in Gaussian % check whether left side of Gaussian extends below xmin: if inx < ng+1 % adjust start indices accordingly: k0 = 1; - g0 = ng-inx+1; + g0 = ng-inx+2; end % check whether right side of Gaussian extends above xmax: if inx > length(kd)-ng @@ -34,9 +34,15 @@ for i = 1:length(data) % for every data value ... % add Gaussian on kernel density: kd(k0:k1) = kd(k0:k1) + kernel(g0:g1); end -kd /= length(data); % normalize by number of data points +kd = kd/length(data); % normalize by number of data points -% plot kernel density: -plot(x, kd) +% plot the computed kernel density: +plot(x, kd, 'b', 'linewidth', 4, 'displayname', 'manual') +hold on +% use the ksdensity() function instead: +[kd, x] = ksdensity(data, x, 'Bandwidth', sigma); +plot(x, kd, '--r', 'linewidth', 4, 'displayname', 'ksdensity()') +hold off xlabel('x') ylabel('Probability density') +legend('show') \ No newline at end of file diff --git a/statistics/lecture/kerneldensity.py b/statistics/lecture/kerneldensity.py index 6f93820..cfdfb62 100644 --- a/statistics/lecture/kerneldensity.py +++ b/statistics/lecture/kerneldensity.py @@ -11,7 +11,7 @@ def kerneldensity(data, xmin, xmax, sigma=1.0) : dx = 0.05*sigma xg = np.arange(-4.0*sigma, 4.0*sigma + 0.5*dx, dx) gauss = np.exp(-0.5*xg*xg/sigma/sigma)/np.sqrt(2.0*np.pi)/sigma - ng = len(gauss)/2 + ng = len(gauss)//2 x = np.arange(xmin, xmax+0.5*dx, dx) kd = np.zeros(len(x)) for xd in data: