%
%  math495/stat490 -- 24 nov 03 -- djm
%
%  w13rw.m:  

clear;

%  walk parameters
T = 1.0;  sig2 = 1.0;
Nst = 300;  dt = T/Nst;  dx = sqrt(sig2*dt);  Nwalks = 10000;

%  initialize random number generator
phone = 2914814;
rand('state',phone);  %  initialize with my office phone number

%  10 sample random walks
figure(1);  clf

walks = cumsum(2*(rand(Nst,10)>0.5)-1)*dx;
plot(dt*(1:Nst),walks)

title('\bf 10 sample random walks')   
ylabel('\bf W(T)-axis')
xlabel('\bf T-axis')

%  ensemble statistics
clear walks
walks = dx*sum(2*(rand(Nst,Nwalks)>0.5)-1);

bins = -4:2*dx:4;
count = hist(walks,bins);

figure(2); clf
plot(bins,count/Nwalks/(2*dx),'*')
hold on
norm = exp(-bins.^2/(2*sig2*T))/sqrt(2*pi*sig2*T);
plot(bins,norm,'r')
title(['\bf probability density of ' num2str(Nwalks) ' rw''s'])   
xlabel('\bf W-axis')
ylabel('\bf f(W)')