This repository has been archived on 2021-05-17. You can view files and clone it, but cannot push or open issues or pull requests.
scientificComputing/statistics/exercises/normprobs.m

58 lines
1.6 KiB
Matlab

%% (a) generate normal distributed random numbers:
n = 10000;
x = randn(n, 1);
%% (b) plot histogram and compare with pdf:
bw = 0.5;
bins=[-5:bw:5];
n = hist(x, bins);
p = n/sum(n)/bw;
subplot(1, 2, 2);
bar(bins, p);
hold on;
xx = [bins(1):0.01:bins(end)];
gauss = exp(-0.5*xx.^2.0)/sqrt(2*pi);
plot(xx, gauss, 'r', 'linewidth', 2);
hold off;
xlim([-5 5])
xlabel('x');
ylabel('p(x)');
%% (c) fraction of data in +/- 1 sigma:
nsig = sum((x>=-1.0)&(x<=1.0));
Psig = nsig/length(x);
fprintf('%d of %d data elements, i.e. %.2f%% are contained in the interval -1 to +1\n\n', ...
nsig, length(x), 100.0*Psig );
%% (d) intgegral over pdf:
dx = 0.01;
xx = -10:dx:10; % x-values
pg = exp(-0.5*xx.^2)/sqrt(2*pi); % y-values Gaussian pdf
% integral over the whole range:
P = sum(pg)*dx;
fprintf( 'Integral over the Gaussian pdf is %.3f\n', P );
% integral from -1 to 1:
P = sum(pg((xx>=-1.0)&(xx<=1.0)))*dx; % we need to use xx, not the random numbers x!
fprintf( 'Integral over the Gaussian pdf from -1 to 1 is %.4f\n\n', P );
%% (e)
for sigma = [1.0, 2.0, 3.0]
Pdata = sum((x>=-sigma)&(x<=sigma))/length(x);
Ppdf = sum(pg((xx>=-sigma)&(xx<=sigma)))*dx;
fprintf( 'Integral over the Gaussian pdf from -%.0f to +%.0f is %.4f\n', sigma, sigma, Ppdf );
fprintf( 'Probability of data in range from -%.0f to +%.0f is %.4f\n\n', sigma, sigma, Pdata );
end
%% (f)
for P = [0.5, 0.9, 0.95, 0.99]
for upper = xx(xx>0.0)
Ppdf = sum(pg((xx>=-upper)&(xx<=upper)))*dx;
if Ppdf > P
fprintf('-%.2f < x < +%.2f: P=%.2f\n', upper, upper, P);
break
end
end
end