close all
clear all
clc

Nh = input('Number of histories per trial: ');

width = 1/sqrt(Nh);
disp(sprintf('Gaussian width = %g',width))

Nt = input('Number of trials: ');

t = zeros(1,Nt);
for i = 1:Nt
   mean = 0;
   for n = 1:Nh
      if (rand < 1/2)
         mean = mean - 1;
      else
         mean = mean + 1;
      end
   end
   t(i) = mean/Nh;
end

tbar = sum(t)/Nt;
s2t  = sum((t - tbar).^2)/(Nt -1); % Width of the distribution
st = sqrt(s2t);
s2tbar = s2t/Nt;
disp(sprintf('Estimated width of the distribution = %g',st))
disp(sprintf('Estimated mean = %g +/- %g ',tbar,sqrt(s2tbar)))

M = 64; % Mesh of display grid
T = zeros(1,M);
tmin = min(t);
tmax = max(t);
x = linspace(tmin, tmax, M);
a = (tmax-tmin)/(M-1);
b = tmin - a;

for i = 1:Nt
   index = ceil((t(i) - b)/a);
   if (index < 1 || index > M)
      disp(sprintf('Index = %g is out of range',index))
   else
      T(index) = T(index) + 1;
   end
end
T = M*T/Nt;
% Integrate[Exp[-(x/s)^2],{x,-Infinity,+Infinity}]
plot(x,T,'ko')