teaching/scientificComputing
Archived
3
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

some code fixes

Browse Source
This commit is contained in:
Jan Benda 2017-11-27 10:26:03 +01:00
parent 82f40834ea
commit 4e17882dd5
3 changed files with 28 additions and 16 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
statistics/code/cumulative.out Normal file
View File

@ -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

36
statistics/code/gaussiankerneldensity.m
View File

@ -1,15 +1,15 @@
data = randn(100, 1); % generate some data data = randn(100, 1); % generate some data
sigma = 0.2; % standard deviation of Gaussian kernel sigma = 0.2; % std. dev. of Gaussian kernel
xmin = -4.0; % minimum x value for kernel density xmin = -4.0; % minimum x value for kernel density
xmax = 4.0; % maximum x value for kernel density xmax = 4.0; % maximum x value for kernel density
dx = 0.05*sigma; % step size 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 xg = [-4.0*sigma:dx:4.0*sigma]; % x-axis for single Gaussian kernel
% single Gaussian kernel: % single Gaussian kernel:
kernel = exp(-0.5*(xg/sigma).^2)/sqrt(2.0*pi)/sigma; kernel = exp(-0.5*(xg/sigma).^2)/sqrt(2.0*pi)/sigma;
ng = (length(kernel)-1)/2; % half the length of the Gaussian ng = floor((length(kernel)-1)/2); % half the length of the Gaussian
x = [xmin:dx:xmax+0.5*dx]; % x-axis for kernel density x = [xmin:dx:xmax+0.5*dx]; % x-axis for kernel density
kd = zeros(1, length(x)); % vector for kernel density kd = zeros(1, length(x)); % vector for kernel density
for i = 1:length(data) % for every data value ... for i = 1:length(data) % for every data value ...
xd = data(i); xd = data(i);
% index of data value in kernel density vector: % index of data value in kernel density vector:
inx = round((xd-xmin)/dx)+1; inx = round((xd-xmin)/dx)+1;
@ -17,13 +17,13 @@ for i = 1:length(data) % for every data value ...
k0 = inx-ng; k0 = inx-ng;
% end index for Gaussian in kernel density vector: % end index for Gaussian in kernel density vector:
k1 = inx+ng; k1 = inx+ng;
g0 = 1; % start index in Gaussian g0 = 1; % start index in Gaussian
g1 = length(kernel); % end index in Gaussian g1 = length(kernel); % end index in Gaussian
% check whether left side of Gaussian extends below xmin: % check whether left side of Gaussian extends below xmin:
if inx < ng+1 if inx < ng+1
% adjust start indices accordingly: % adjust start indices accordingly:
k0 = 1; k0 = 1;
g0 = ng-inx+1; g0 = ng-inx+2;
end end
% check whether right side of Gaussian extends above xmax: % check whether right side of Gaussian extends above xmax:
if inx > length(kd)-ng if inx > length(kd)-ng
@ -34,9 +34,15 @@ for i = 1:length(data) % for every data value ...
% add Gaussian on kernel density: % add Gaussian on kernel density:
kd(k0:k1) = kd(k0:k1) + kernel(g0:g1); kd(k0:k1) = kd(k0:k1) + kernel(g0:g1);
end end
kd /= length(data); % normalize by number of data points kd = kd/length(data); % normalize by number of data points
% plot kernel density: % plot the computed kernel density:
plot(x, kd) 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') xlabel('x')
ylabel('Probability density') ylabel('Probability density')
legend('show')

2
statistics/lecture/kerneldensity.py
View File

@ -11,7 +11,7 @@ def kerneldensity(data, xmin, xmax, sigma=1.0) :
dx = 0.05*sigma dx = 0.05*sigma
xg = np.arange(-4.0*sigma, 4.0*sigma + 0.5*dx, dx) 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 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) x = np.arange(xmin, xmax+0.5*dx, dx)
kd = np.zeros(len(x)) kd = np.zeros(len(x))
for xd in data: for xd in data: