function randomwalk(N,p, M)


x=zeros(1,M);

for i=1:M
  xi = rand(1,N);
  
  Np = length(xi(xi<p));
  x(i) = 2*Np - N;

  %for j=1:N
  %  if (xi(j)<p) xi(j) = +1; else xi(j) = -1; end
  %end
  %x(i) = sum(xi);
  
end

xmean = mean(x)
xvar = std(x)^2

bin = -N:2:N;
h = hist(x,bin);
Pexp = h/M/2;

figure(1)
x=-100:0.01:100;

Np = (N+x)/2;
Nm = (N-x)/2;

P = p.^Np .* (1-p).^Nm * gamma(N+1) ./ (gamma(Np+1) .* gamma(Nm+1))/2;
plot(x,P,'b');
hold on
stairs(bin-1,Pexp,'r');
hold off

xmean_thr = N*(2*p-1);
xvar_thr = 4*N*p*(1-p);

xmean_thr = xmean;
xvar_thr = xvar;

Pgauss = 1/sqrt(2*pi*xvar_thr)*exp(-0.5*(x-xmean_thr).^2/xvar_thr);

hold on
plot(x,Pgauss,'m:')
hold off
