樓上兩位的回答都很用心，也很精彩，贊壹個。
我的代碼主要有以下優點：
（1）用稀疏矩陣存儲a，克服內存不足問題（N取100萬，使用的內存還不到20M）。
（2）繪圖動態顯示N次模擬過程中r/R的變化。
代碼如下（同時已作為附件上傳）：
N?=?1000000;M?=?2*N;
a?=?sparse(M+1,?M+1);
j?=?N?+?1;
k?=?N?+?1;
b?=?ceil(4*rand(1,N));
%?繪圖顯示計算過程(為提高效率，每n次循環輸出壹個點)
n?=?500;
v?=?zeros(1,?fix(N/n)+1)?*?NaN;
h?=?plot(1?:?fix(N/n)+1,?v,?'b-',?NaN,?NaN,?'ro');
xlabel('N');
ylabel('r/R');
set(gcf,?'DoubleBuffer',?'on');
set(gca,?'Xlim',?[1?fix(N/n)+1]);
t=ceil(exp(1:15));
set(gca,'xtick',t/n,'xgrid','on','xticklabel',t)
for?i?=?1:N
switch?b(i)
case?1,?j=j+1;?a(j,k)=a(j,k)+1;
case?2,?k=k+1;?a(j,k)=a(j,k)+1;
case?3,?j=j-1;?a(j,k)=a(j,k)+1;
case?4,?k=k-1;?a(j,k)=a(j,k)+1;
end
%?更新繪圖
if?~rem(i,?n)?||?i?==?N
%?註意：不能用?sum(a(:)~=0)?進行計算，否則容易導致內存不足
R?=?full(sum(sum(a~=0)));
c?=?fix(log(i));
r?=?full(sum(sum(a==c)));
idx?=?fix(i/n)?+?1;
v(idx)?=?r/R;
set(h(1),?'yData',?v);
set(h(2),?'xData',?idx,?'yData',?v(idx));
title(['N?=?'?int2str(i)]);
drawnow
end
end
%?輸出結果
fprintf('R=%i,?r=%i,?r/R=%.3g\n\n',?R,?r,?r/R);
%?統計各方向移動的次數（驗證隨機數的均勻性）
for?i?=?1?:?4
fprintf('方向%i的次數為%i\n',?i,?sum(b==i));
end
某次程序的輸出如下：
R=204146,?r=2856,?r/R=0.014方向1的次數為249327
方向2的次數為250179
方向3的次數為250085
方向4的次數為250409
簡單說明幾點：
1、由於是隨機模擬，每次運行的結果都會有差別。
2、移動是壹個前後關聯的過程，所以隨機數序列不僅要求均勻，還應該獨立（相鄰的隨機數之間不相關）。
3、圖中的虛線表示 int(ln(n)) 發生變化的N，計算格點次數變了，所以通常表現為不連續。
4、從圖中的變化趨勢看，沒有收斂到某壹個數的明顯跡象，我運行兩次的結果分別是0.0118和0.014。
5、在我的機器上，取N=100萬，運行壹次的時間大約是10分鐘。
6、程序對MATLAB版本沒有特別要求，在6.5、2008a、2012b上測試過，都可以正常運行。