/* Program to accompany the article "The dunkard's walk in 2-D," Rick Wicklin, published 12AUG2015, The DO Loop blog, http://blogs.sas.com/content/iml/2015/08/05/2d-drunkards-walk.html */ proc iml; call randseed(54321); dim = 2; /* number of directions */ NumSteps = 10; /* number of steps */ coord = j(NumSteps,1,.); /* allocate vector for random direction */ call randgen(coord, "Table", j(1,dim,1)/dim); /* random directions */ incr = sample({-1 1}, 1 // NumSteps); /* random step +/-1 */ step = T(1:NumSteps); print step coord incr; x = j(NumSteps, dim, 0); /* (N x d) array */ idx = sub2ndx(dimension(x), step||coord); x[ idx ] = incr; print x; /***************************************************/ proc iml; start vcusum(x); c = j(nrow(x), ncol(x), .); do i = 1 to ncol(x); c[,i] = cusum(x[,i]); end; return (c); finish; start DrunkardWalk(dim, NumSteps); coord = j(NumSteps,1,.); /* allocate vector for random direction */ call randgen(coord, "Table", j(1,dim,1)/dim); /* random directions */ incr = sample({-1 1}, 1 // NumSteps); /* random step +/-1 */ step = T(1:NumSteps); x = j(NumSteps, dim, 0); /* (N x d) array */ idx = sub2ndx(dimension(x), step||coord); x[ idx ] = incr; return( x ); finish; store module=_all_; call randseed(54321); dim = 2; /* number of directions */ NumSteps = 10; /* number of steps */ x = DrunkardWalk(dim, NumSteps); Y = vcusum(j(1,dim,0) // x); /* show initial position at t=0 */ Q = t(0:nrow(x)) || Y; create walk from Q[c={"Step" "x" "y"}]; append from Q; close; quit; title "A Single Drunkard's Walk"; proc sgplot data=walk; series x=x y=y / datalabel=Step markers; refline 0/ axis=x; refline 0/ axis=y; xaxis integer min=-2 max=2; yaxis integer min=-2 max=2; run; /************************************************/ proc iml; load module=(vcusum DrunkardWalk); dim = 2; NumSteps = 6; call randseed(1234); NumWalks = 10000; finalPos = j(NumWalks, dim); do i = 1 to NumWalks; x = DrunkardWalk(dim, NumSteps); finalPos[i,] = x[+,]; end; /* summarize and visualize 2-D data */ create walk from finalPos[c={"x" "y"}]; append from finalPos; close; submit; proc freq data=walk noprint; /* compute frequency distribution */ tables x*y / list sparse out=freqOut; run; endsubmit; /* read into IML vectors */ use freqOut; read all var {x y count}; close; nX = range(x) + 1; nY = range(y) + 1; W = shape(count, nX, nY); /* reshape into rectangular matrix */ ods graphics / width=600px height=500px; title "10,000 Drunkard's Walks"; title2 "Distribution of Counts for Final Positions"; call HeatmapCont(W) xValues = unique(x) yValues=unique(y); pct = W / sum(W); print pct; quit;