ntime = 100;
npart = 1000;
steps=2*randi(2,ntime,npart)-3;
P=cumsum(steps);
figure(1)
plot(P)
ylabel('Position')
xlabel('Step Count')
title('1-D Random Walk')
figure(6), subplot(2,1,1), plot(steps(:,1),'o-')
figure(6), subplot(2,1,2), plot(P(:,1),'o-')
%%
% Histogram at ntime
figure(2)
h=histogram(P(ntime,:),10);
xlabel('Distance, steps')
ylabel('Counts')
ydist=h.Values/npart/h.BinWidth;
xdist=h.BinEdges(1:end-1)+h.BinWidth/2;
%%
figure(3)
plot(xdist,ydist,'.-k')
xlabel('Distance, steps')
ylabel('Probability')
title(['Probability distribution after ',num2str(ntime),' steps'])
%%
% Gaussian distribution of width 
x=-ntime:1:ntime;
mu=0;
sigma=std(P(ntime,:));
eta=(x-mu)/sigma/sqrt(2);
G=exp(-eta.^2)/(sigma*sqrt(2*pi));
figure(3)
hold on
plot(x,G,'-r')
axis([-ntime/2 ntime/2 0 max(G)])
hold off
%%
figure(4)
Ps=sort(P(ntime,:));
plot(Ps,(1:npart)/npart,'-k')
xlabel('Distance, X')
ylabel('Cumulative probability P(x>X)')
title('Cumulative probability distribution after 100 steps')
%%
% Integral of Gaussian distribution of width 10
E=(1+erf(eta))/2;
figure(4)
hold on
plot(x,E,'-r')
axis([-ntime/2 ntime/2 0 1])
hold off
%%
% Mean square displacement
SD=P.^2;
MSD=mean(SD,2);
figure(5)
t=1:ntime;
[ b,stdb ]=linregr(t',MSD);
y=t*b;
plot(t,MSD,'ok',t,y,'-r')
text(ntime/2,max(MSD)/2,['slope = ',num2str(b),'+/-',num2str(stdb)])