Matlab Monte Carlo

Buffon’s needle problem

pi ~ (2l/a)(m/n)

function buffonNeedle(a,l,q)
format long;
n=0;
h=rand(q,1);
h=(a/2).*h;
theta=rand(q,1);
theta=(pi).*theta;
for m=1:q
    if h(m)<=sin(theta(m))*(l/2)
        n=n+1;
    end
end
(2*l*m)/(a*n)
%save as buffonNeedle.m Note that
%a is leading; l is length of needle; q is repeated times;

When executing, gives

buffonNeedle(10,6,1e6)

ans =

   3.141385927114611;

答案最大误差约为0.1，或者千分之5以内。

Area of the unit circle

function ucpi(n)
format long;
x=rand(n,1);
y=rand(n,1);
i=1;
m=0;
for i=1:n
    if x(i)^2+y(i)^2<=1
        m=m+1;
    end
end
4*m/n

%save as ucpi.m

Gives:

ucpi(1e6)

ans =

   3.14282400000000

写这个程序的时候，我发现一个很有意思的事情。一开始，我是通过下面这种方法实现的：

%THIS FUNCTION DOES NOT WORK AS INTENDED----------------------------------

function ucpi(n)
format long;
t=rand(n,1);
i=1;
m=0;
for i=1:n
    if t(i)<=sqrt(1-t(i)^2)
        m=m+1;
    end
end
4*m/n

很好理解我是怎么想的，选取一个数，比如0.2，带到单位圆里算，如果他在单位圆里面的话，就会满足x^2 + y^2 <1；
或者， y<= sqrt(1-x^2) 也就是t(i)<=sqrt(1-t(i)^2)这一行。
执行，发现答案在2.8附近，也就是正确答案pi的0.9倍。
怎么回事？
整理这一行公式：(2t(i)^2<=1)时计数，和源程序的区别，在于此时x和y兼并成了一个t。
也就是说，这个程序所给出的 m/n的数值，为x=y这条线上，在圆内和在圆外部分的占比，或者， 1/sqrt(2)；
那试着跑一下：

function ucpi(n)
format long;
t=rand(n,1);
i=1;
m=0;
for i=1:n
    if t(i)<=sqrt(1-t(i)^2)
        m=m+1;
    end
end
m/n

>>>>>>>--------------------------------------------------------------------
ucpi(1e6)

ans =

   0.70775100000000
%     0.70710678118655 为 1/sqrt(2)的值

哦，原来如此。哈哈哈。