(* simple probability distribution *)
rho[x_]:=If[x<0||x>1,0,1];
(* empty array storing the markov chain *)
xlist=ConstantArray[0,ntrials];
(* start x at 0.5, do 10^6 trials, and  *)
xcur=0.5;
ntrials=10^6;
delta=0.05;
propose:=RandomReal[{-delta,delta}];

(* populate the Markov chain *)
Do[
xprop=xcur+propose;
If[RandomReal[]<rho[xprop]/rho[xcur],xcur=xprop];
xlist[[i]]=xcur; (* Could also write AppendTo[xlist,xprop] but it would be slower. *)
,{i,1,ntrials}];

ListLinePlot[xlist[[1;;1000]],PlotLabel->"Part of the markov chain",PlotRange->{-0.01,1.01}]

(* binning error calculation. Note the extra factor of M \[Equal] Length[bins] *)
binnedError[binsize_]:=
Module[{bins,abar},
bins=Mean/@Partition[xlist,binsize];
abar=Mean[bins];
Sqrt[Total[(bins-abar)^2]/Length[bins]^2]
];

(* Plot the error estimates over increasing bin size. *)
errorlist=Table[{i,binnedError[i]},{i,1,10000,10}];
maxerror=Max[errorlist[[All,2]]];
Show[
ListLinePlot[errorlist,PlotRange->All,AxesLabel->{"bin size","error"},
PlotLabel->"Binned Error Estimates of Mean[xlist]="<>ToString[Mean[xlist]]],
Plot[maxerror,{i,1,10000},PlotStyle->{Thick,Dashed,Red}]
]

(* Plot the autocorrelation and estimate the autocorrelation time *)
autocorrelation[i_]:=Mean[xlist[[1;;-i]]xlist[[i;;]]]-Mean[xlist[[1;;-i]]]Mean[xlist[[i;;]]];
autos=Table[{i,autocorrelation[i]},{i,1,2000,10}];
(* fit it to an exponentially decreasing curve. We discard c but we need it for the fit *)
fit=FindFit[autos,c E^(-x/tau),{c,tau},x];
(* Plot everything with fancy formatting. *)
Show[
Plot[c E^(-n/tau)/.fit,{n,0,2000},PlotRange->All,PlotLegends->{"fit"},
PlotStyle->Directive[Dashed,Thin,Red],AxesLabel->{"n","correlator"}],
ListPlot[autos,PlotRange->All,PlotLegends->{"data"}],
PlotLabel->"Autocorrelation time estimated to be "<>ToString[tau/.fit]
]