ビュフォンの針のシミュレーション

格子に針を落とした時の針と講師が交わる確率を求めようとして下のコードを書きましたが、理論値と大きくずれてしまいます。

MATLAB初心者で見様見真似で作りました。コードのどこが間違っているか教えていただきたいです。(短いコードなのですべて投稿させていただきました。)

clear
clc
N = 1000;%針の数
L = 0.20;%針の長さ
xb = L + rand(1,N)*(1-2*L);%針の始点のx座標
yb = L + rand(1,N)*(1-2*L); %針の始点のx座標
angs = rand(1,N)*360;%針の傾き
xe = xb + L*cosd(angs);
ye = yb + L*sind(angs);
ax = axes;
plot(ax,[xb;xe],[yb;ye]) %針の表示
axis square
hold on
glines = 0:L:1;
for i = 1:length(glines) %x軸に平行な線の表示
   xline(ax,glines(i));
end
for p = 1:length(glines) %y軸に平行な線の表示
    yline(ax,glines(p));
end

r = floor(yb/L) ~= floor(ye/L); %y軸に針が交わるか交わらないかの判定(結果は行ベクトル形式で出力)
q = floor(xb/L) ~= floor (xe/L); %x軸に針が交わるか交わらないかの判定(結果は行ベクトル形式で出力)
qr = [r;q]; %2つのベクトルを結合(2列のベクトルにする)
R = sum(r);%y軸に交わった針の個数
Q = sum(q);%x軸に交わった針の個数
zeroqr = all(qr,1);
A = sum(zeroqr);%x軸にとy軸共に交わった針の個数
probability = (R+Q-A)/N%針が講師の交わる確率

probability = 0.9550

Theoreticalvalue =1-(1-2/pi)^2 %理論値を出力

Theoreticalvalue = 0.8680

[Atsushi Ueno]

ビュフォンの針 - Wikipedia：縦線と交わる確率、横線と交わる確率はそれぞれ で理論値通りです。

Theoreticalvalue =1-(1-2/pi)^2 %理論値を出力

自分も始めは表の様に考えました。プログラム中の上記理論値は黄色で塗りつぶした「針と(縦横片方/両方の)格子線が交わる確率」ですよね。但し、これは縦横各々の発生確率が互いに独立している場合の考え方です。交差確率には針の角度が関係していて「縦の格子と交わり難い場合は横の格子と交わり易い」事になるので、相互作用を考慮した確率分布を考える必要があります。

...どうすれば良いのかな🤔

[Atsushi Ueno]

Estimating pi Using Buffon’s Method 著者 Loren Shure, August 22, 2019