%% (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