在一個圓 A 內隨機選三個點，這三個點的外接圓 B 包含於圓 A 的概率是多少？

09-02

就是隨機產生點，只取圓A內的點為樣本

怎麼隨機產生點

--------------------

如果這是個題那可以算一算,不過如果是自己沒事想出來的那還是別算了,因為多半都沒有初等表達.

1.用蒙特卡洛跑出來結果在0.4左右，理論方法不會

2.取點的是 $(sqrt{u}cos{ heta},sqrt{u}sin{ heta})$ , $u$ 和 $heta$ 都是均勻分布。用其他方法結果也會不同

（考慮到這個問題與Mathematica無關，把結論放在最前，下面的內容對MMA不感興趣的同學們可以無視了）

用Mathematica隨便寫了個蒙特卡洛，幾乎全部代碼來自於內置幫助

In[141]:= UniformDisk /: Random`DistributionVector[UniformDisk[r_?Positive], n_Integer, prec_?Positive] := r*Sqrt[RandomReal[1, n, WorkingPrecision -&> prec]]* Transpose[{Cos[#], Sin[#]} [ RandomReal[2 Pi, n, WorkingPrecision -&> prec]]]; CircumCenter[{x1_, y1_}, {x2_, y2_}, {x3_, y3_}] := Module[{a, bx, by}, a = Det[{{x1, y1, 1}, {x2, y2, 1}, {x3, y3, 1}}]; bx = -Det[{{x1^2 + y1^2, y1, 1}, {x2^2 + y2^2, y2, 1}, {x3^2 + y3^2, y3, 1}}]; by = Det[{{x1^2 + y1^2, x1, 1}, {x2^2 + y2^2, x2, 1}, {x3^2 + y3^2, x3, 1}}]; {-bx/(2 a), -by/(2 a)}]; With[{n = 5*10^4}, Sum[Module[{p1, p2, p3, c}, {p1, p2, p3} = RandomReal[UniformDisk[1], 3]; c = CircumCenter[p1, p2, p3]; Boole[Norm[c] + Norm[p1 - c] &<= 1] ], {n}]/n // N] // AbsoluteTiming Out[143]= {6.723612, 0.40096}

第一段是定義在圓盤上取隨機點，來自於幫助的tutorial/RandomNumberGeneration#132002727部分

第二段是定義計算三角形外接圓圓心的函數，來自於Circle的幫助中「應用」-&>「定義三角形外接圓中心」部分

第三段是計算蒙特卡洛的部分，這裡取5*10^4個三角形計算

最後運行將近7秒，得到結果0.4左右

這個代碼最大的優點就是簡單無腦，幾乎全是從幫助里複製粘貼的；最大的缺點就是慢，畢竟幫助里的例子不是為了大規模計算設計的沒考慮過效率優化。

要提速的話，大概可以從向量化、編譯、並行三方面下手。我稍微試了下並行+編譯能提速10倍左右，但代碼變得丑且巨長，就不發上來了。向量化嫌麻煩就沒做，哪天有功夫回來研究下

雖然我覺得這種純簡單粗暴的程序再怎麼優化也比不上直接用C。。。。。。

======================15年9月17日更新=============================

回答 @朱裕德的問題，字數超了沒法貼到評論里，就直接更新答案了

首先測試你下你的代碼的速度

In[1]:= n = 10000; Count[Table[ circle = Circumsphere[ RandomPoint[Disk[], 3]]; (Norm[circle[[1]]] + circle[[2]]) &<= 1, n], True]/n // N // AbsoluteTiming Out[2]= {2.06577, 0.4001}

Circumsphere是內置函數不支持編譯，想通過編譯提速的話需要手寫一個求三點外接圓的函數，雖然麻煩點但是確實能提速不少。但是你這個慢的主要原因其實不在Circumsphere而在RandomPoint：

In[2]:= pt = RandomPoint[Disk[], {10^4, 3}]; // AbsoluteTiming


Out[2]= {0.00132639, Null}
In[3]:= Do[RandomPoint[Disk[], 3], {10^4}]; // AbsoluteTiming

Out[3]= {1.7027, Null}

可以看出一次生成10^4個隨機點比生成1個隨機點10^4次快了1000倍以上。。。。。你的代碼里大部分時間都是在這裡浪費掉的。修改以後

In[58]:= AbsoluteTiming[ ptlist = RandomPoint[Disk[], {10^4, 3}]; Count[(Norm[#[[1]]] + #[[2]]) &<= 1 /@ Circumsphere /@ ptlist, True]/10^4 // N ] Out[58]= {0.162775, 0.4072}

在我的電腦上比你的代碼快了15倍左右

然後就是編譯和向量化的提速了。直接把上面CircumCenter定義里的行列式都展開帶入進去，再稍微整理一下即可得到圓心的表達式

In[4]:= AbsoluteTiming[ cf = Compile[{{ptlist, _Real, 2}}, Module[{x1, x2, x3, y1, y2, y3, center, r}, {{x1, y1}, {x2, y2}, {x3, y3}} = ptlist; center = {(x3^2 (y1 - y2) + x1^2 (y2 - y3) + x2^2 (-y1 + y3) - (y1 - y2) (y2 - y3) (-y1 + y3))/(2 (x3 (y1 - y2) + x1 (y2 - y3) + x2 (-y1 + y3))), ((x1 - x2) (x2 - x3) (-x1 + x3) - (x2 - x3) y1^2 - (-x1 + x3) y2^2 - (x1 - x2) y3^2)/(2 (x3 (y1 - y2) + x1 (y2 - y3) + x2 (-y1 + y3)))}; r = Norm[{x1, y1} - center]; Norm[center] + r &<= 1], RuntimeAttributes -&> {Listable}]; ptlist = RandomPoint[Disk[], {10^4, 3}]; Count[cf@ptlist, True]/10^4 // N]

Out[4]= {0.00633942, 0.4039}

可以看到即使加上編譯花費的時間也比原先快了20倍以上

事實上如果加上CompilationTarget -&> "C"選項，單純的計算速度還能加快將近一倍，只不過相應的編譯所需時間也會有所增加，對於10^4這個數量級就有些得不償失了，大約要10^6以上的量級才能補回編譯所需的時間

如圖所示

最後放張n=10^8時的測試截圖

用時40秒，所得結果比上面的更接近0.4了，所以我有充分的理由懷疑本題可以用理論方法求出解析解，概率就是2/5，但是實在不會證。。。。。。。。

@無影東瓜，實際上通過一些 meta-programming，代碼可以很簡潔。

Mathematica:

f=Evaluate[UnitStep[1-Norm@#-#2]@@Circumsphere@{{#1,#2},{#3,#4},{#5,#6}}]; Total[f@@RandomPoint[Disk[],{1*^6,3}]~Flatten~{{2,3},{1}}]

(* 400164 *)

待續...

寫了這麼個程序，我感覺我這是在抹黑mma

n = 2000; Count[Table[ RegionIntersection[r = Circumsphere[RandomPoint[Disk[], 3]], Disk[]] === Circle @@ r, n], True]/n // N

｛0.401｝

For的效率更是低到丟人

For[i = 0; sum = 0, i &<= 1000, i++, If[RegionIntersection[r = Circumsphere[RandomPoint[Disk[], 3]], Disk[]] === Circle @@ r, sum++]] sum/1000.

｛0.376｝

慢的原因是沒有想到快的方法判斷一個Region是不是屬於另一個Region，只能用交來弄，所以慢死了

按 @秦吉寧的提議改了改，快了大概10倍左右，不過這裡情況特殊，都是圓之間的判斷，所以可以這麼改，但mma不能判斷一個region是不是屬於另一個region的事我已經碰到好幾次了，不容易處理，取交集來判斷慢到掉渣，而且這題這麼改速度問題還是有點著急，編譯提速也沒學好，一編就報錯，待我@幾個編譯大師出來指教一下，代碼如下：

n = 10000; Count[Table[ circle = Circumsphere[ RandomPoint[Disk[], 3]]; (Norm[circle[[1]]] + circle[[2]]) &<= 1, n], True]/n // N

｛0.4008｝

@無影東瓜@蘋果

我好像能給個計算表達式。

這裡假定隨機是圓內關於面積的平均分布。

我們給這三個點一個順序，x1,x2,x3，那麼我們可以旋轉坐標使得x1x2與x軸坐標平行：

那麼我可以把這兩點的坐標設為(x1,t),(x2,t)

那麼現在放置x3，設f(x1,x2,t)為三個點的外接圓 B 在圓 A內的概率。

f(x)如圖上白色部分所示：

其中上下兩個圓與大圓相切

那麼我們所求的就是一個積分式：

$int^{1}_{-1}{frac{2 sqrt{1-t^2}}{pi}int_{-sqrt{1-t^2}}^{sqrt{1-t^2}}{int_{-sqrt{1-t^2}}^{sqrt{1-t^2}}{frac{f(x_1,x_2,t)}{4(1-t^2)}dx_1}dx_2}dt}$

設x1& $int^{1}_{-1}{frac{4 sqrt{1-t^2}}{pi}int_{-sqrt{1-t^2}}^{sqrt{1-t^2}}{int_{x_1}^{sqrt{1-t^2}}{frac{f(x_1,x_2,t)}{4(1-t^2)}dx_2}dx_1}dt}$

其中關於t的項是相應的概率密度（可能有問題，我沒想通）。

然後我們來求這個f(x):

來求這兩個圓的圓心的縱坐標，由於圓與大圓想切，圓心 $(frac{x_1+x_2}{2},s)$ 滿足的方程為：

小圓半徑＋大小圓心間距離＝大圓半徑

$sqrt{frac{x_1-x_2}{2}^2+(t-s)^2}+sqrt{frac{x_1+x_2}{2}^2+s^2}=1$

這裡設中間變數：

$frac{x_2-x_1}{2}=a>0,frac{x_1+x_2}{2}=b$

帶入mathematica大法算一算得：

$ext{Solve}left[sqrt{a^2+(s-t)^2}+sqrt{b^2+s^2}=1,s ight]$

$s_1 o frac{-a^2 t-sqrt{a^4-2 a^2 b^2+2 a^2 t^2-2 a^2+b^4+2 b^2 t^2-2 b^2+t^4-2 t^2+1}+b^2 t+t^3-t}{2 left(t^2-1 ight)}$

$s_2 o frac{-a^2 t+sqrt{a^4-2 a^2 b^2+2 a^2 t^2-2 a^2+b^4+2 b^2 t^2-2 b^2+t^4-2 t^2+1}+b^2 t+t^3-t}{2 left(t^2-1 ight)}$

那麼圓心坐標（b,s1),(b,s2)

然後就可以算弓形面積

$S_{2}=-( a(s_{2}-t) )+(a^2+(t-s_2)^2)2tan^{-1}{frac{a}{s_2-t} }$

$S_{1}=-( a(t-s_{1}) )+(a^2+(t-s_1)^2)2tan^{-1}{frac{a}{t-s_1} }$

$f(x)=frac{S-2(S_1+S_2)}{pi}$

其中S為兩圓面積之和。噗，算出來不對啊，0.211467，那麼說明方法不對，原來的面積等概率不能化為那個積分式子。

我從Springer出版社ULNP（Undergraduate Lecture Notes in Physics）系列中的Inside Interesting Integrals（作者Paul J.Nahin）一書中讀到了以下思路，侵刪。Paul J.Nahin在書中寫道，這一證法由Joseph Edwards在A Treatise on the Integral Calculus一書的817-818頁給出，並且Paul J.Nahin對這一做法有疑問。閱讀時請持保留態度。

在圓C?中隨機且獨立地選取不共線的三點，該三點的外接圓記為C?。C?的半徑記為R，C?的半徑記為r，C?與C?圓心距記為l。

以寬度為dr的圓環包圍C?的周長，使圓周平分該圓環。（描述不嚴謹請包涵。）

圓環的面積

$S=pi (r+frac{dr}{2} )^2-pi (r-frac{dr}{2} )^2=2pi rdr$

則從C?隨機選取一點，該點落在圓環內的概率為

$frac{2pi rdr}{pi R^2}=frac{2r}{R^2}dr$

由於三點的選取互相獨立，故三點都落在圓環內的概率為

$left( frac{2r}{R^2}dr ight) ^3=frac{8r^3}{R^6}(dr)^2dr=frac{8r^3}{R^6}dSdr$

取定積分

$int_{0}^{R-l}int_{0}^{pi R^2}frac{8r^3}{R^6}dSdr$

$=frac{8}{R^6}int_{0}^{2pi }int_{0}^{R-l}int_{0}^{R}r^3ldldrd heta$

$=frac{16pi }{R^6}int_{0}^{R-l}int_{0}^{R}r^3ldldr$

$=frac{4pi }{R^6}int_{0}^{R}(R-l)^4ldl$

做一步換元

$u=R-l$

從而原式

$=frac{4pi }{R^6} int_{0}^{R}u^4(R-u)du$

$=frac{2pi }{15}$

下面引用Paul J.Nahin對該證法的疑問：

How do we really know that one of our admittedly casual manipulations along the way didn』t have a hidden ?aw in it?(Like the claim, for example, that {Δx}2 is the differential area in rectangular coordinates, in which I』ve replaced the Δy in the usual dA=ΔxΔy with another Δx. Or, what about the claim a randomly selected point from the interior of C? is in the Δx band around C?, as that assumes C? is totally inside C??)
我們如何知道在這個過程中，我們那些無疑很隨意的操作沒有隱藏的漏洞呢？比如聲明{Δx}2是平面直角坐標系內面積的微分，而在這一操作中我用Δx替代了常規的面積微分dA=ΔxΔy 中的Δy。再比如，「在C?中隨機選取一點，該點落在包圍C?周長的圓環內」就意味著「C?內含於C?」的聲明又是否正確呢？

在Inside Interesting Integral一書中還給出了MATLAB模擬的結果，代碼如下（保留了注釋，侵刪）：

%Created by P.J.Nahin(9/12/2012)for Inside Interesting Integrals. %This code performs simulations of the circle in a circle problem %one million times. inside=0; for loop1=1:1000000 for loop2=1:3 done=0; while done==0 x=-1+2*rand;y=-1+2*rand; if x^2+y^2&<=1; px(loop2)=x;py(loop2)=y;done=1; end end end ma=(py(2)-py(1))/(px(2)-px(1)); mb=(py(3)-py(2))/(px(3)-px(2)); xc=ma*mb*(py(1)-py(3))+mb*(px(1)+px(2))-ma*(px(2)+px(3)); xc=xc/(2*(mb-ma)); yc=(-1/ma)*(xc-((px(1)+px(2))/2)); yc=yc+(py(1)+py(2))/2; center=sqrt(xc^2+yc^2); radius=sqrt((xc-px(1))^2+(yc-py(1))^2); if center+radius&<=1 inside=inside+1; end end inside/loop1

按照MATLAB模擬實驗一百萬次的結果，頻率大致在0.3990到0.4007(親測)間波動，感覺實驗結果強烈暗示正確的概率是0.4（前面幾位答主也提到這個答案了），而上述解法的理論值2π/15≈0.418879，誤差大概在4.3%~4.7%，感覺不太靠譜。

大概就這樣，如果這個解法確實錯了，歡迎大家在評論區指錯，我會及時更正的_(:3」∠)_

這個題，用蒙特卡洛方法解答的話，我是這樣考慮的：

三個點P1(x1,y1),P2(x2,y2),P3(x3,y3)可以看做一個6維線性空間O中均勻分布的一個點，三點外接圓C2內含於指定圓C1是一個約束，這個約束在6維空間中的體積是V1，而三點內含於C1也是個約束，這個約束在O中的體積為V2，所求概率P=(V1∩V2)/V2。

V2好說，但是V1∩V2的交界應該比較難處理，因為交界是個五維流形，還要一層層的切下去，雖然理論上只要細心能夠一層層的積出來，但是我相信應該有更簡單的辦法。

等等我錯了，V1是V2的真子集，P=V1/V2 ，沒有交界的事情了，直接積V1的體積就好，不過這個約束形式很繁瑣啊

我的解析幾何就是渣，不知道給定坐標的三個點如何計算其外接圓的圓心坐標和半徑。

剩下的事情不難辦，蒙特卡洛法吧。

均勻分布叫隨機，高斯分布也是隨機……隨機的方式太多啦

**********************

如果是等概率選點，直覺告訴我，概率為 0