Mathematica/均勻球面分佈

在本教程中，我們將繪製一組隨機的、均勻分佈的球體表面上的點。這看起來像一個微不足道的任務，但我們會看到“明顯”的解決方案實際上是錯誤的。我們將從這種錯誤方法開始，然後改進它使其正確。

錯誤方法

我們使用的球座標系。

為了在球體表面上繪製一個點，我們需要兩個數字：餘緯度，φ 和經度，θ。我們可以使用緯度，它由 90°−φ 給出，但餘緯度是處理的傳統量。

現在，我們將做一個簡單的假設，我們可以對餘緯度和經度使用均勻分佈。我們只會在半個球體上繪製 6000 個點，因為這將使我們能夠更清楚地看到結果。

讓我們定義這些分佈

theta := RandomReal[ {0, Pi} ]; (* Uniform on [0,Pi], longitude, only goes halfway round*)
phi   := RandomReal[ {0, Pi} ]; (* Uniform on [0,Pi], colatitude, from pole to pole *)

您可以透過更改覆蓋的角度輕鬆更改{0, Pi}到您想要的邊界。

請注意，我們使用了“延遲設定”運算子 (:=）。我們需要這個，因為當我們在下一步構建一個值的表時，我們希望Random[]函式在每個條目中被呼叫。如果我們不這樣做，的值theta和phi將在上面設定並永遠不會改變，導致 6000 個相同的點。我們現在製作資料表

dataPoints = Table[
   {phi, theta},
   {6000}];

讓我們看一下點的分佈

rectPlot = ListPlot[
  dataPoints,
  AxesLabel -> {"\[Phi]", "\[Theta]"},
  Ticks -> {
    {0, Pi/4, Pi/2, 3 Pi/4, Pi},
    {0, Pi/4, Pi/2, 3 Pi/4, Pi}}
  ]

這將產生以下輸出。我們看到，這些點確實似乎是均勻分佈的。

現在，我們將這些球座標對映到笛卡爾空間。這種變換是標準的

$x=r\sin \phi \cos \theta \,$

$y=r\sin \phi \sin \theta \,$

$z=r\cos \phi \,$

為了我們的目的，我們將設定 r=1。我們製作一個名為sphericalToCartesian的函式，用於轉換每個雙元素座標{φ,θ}在球面系統中，到三元素座標，{x,y,z}，在笛卡爾系統中。

sphericalToCartesian[{phi_, theta_}] :=
  {Sin[phi]*Cos[theta],
   Sin[phi]*Sin[theta],
   Cos[phi]};

我們有一個對映，所以剩下的就是使用“對映”運算子來應用它，/@，它將在每個點上執行此操作dataPoints:

mappedPoints = sphericalToCartesian /@ dataPoints;

現在我們有了笛卡爾系統中的一組點，我們想看看它們。我們製作了一組位於每個點的點mappedPoints並顯示它。這裡，我們有預設檢視，以及來自頂部和正面的正投影檢視

 points = Graphics3D[Point /@ mappedPoints];

 Show[points]
 Show[points, ViewPoint->{0,Infinity,0}]
 Show[points, ViewPoint->{0,0,Infinity}]

我們看到，這些點肯定不是均勻分佈的 - 它們在兩極處密度更高。這是因為從球面到笛卡爾座標的對映不保留面積 - 球面空間在對映中被擠壓和壓縮到兩極。為了解決這個問題，我們必須改變我們最初的分佈。

正確方法

為了找到的正確分佈theta和phi，我們必須考慮雅可比行列式對映sphericalToCartesian。這將告訴我們球面空間在每個點上是如何被變換擠壓的。

我們必須首先找到雅可比矩陣。我們暫時保留 r，因為 Mathematica 將對數字進行運算。評估的雅可比矩陣如下所示，但我們將讓 Mathematica 從頭開始計算出來。

$J_{F}(r,\phi ,\theta )={\begin{bmatrix}{\dfrac {\partial x}{\partial r}}&{\dfrac {\partial x}{\partial \phi }}&{\dfrac {\partial x}{\partial \theta }}\\[3pt]{\dfrac {\partial y}{\partial r}}&{\dfrac {\partial y}{\partial \phi }}&{\dfrac {\partial y}{\partial \theta }}\\[3pt]{\dfrac {\partial z}{\partial r}}&{\dfrac {\partial z}{\partial \phi }}&{\dfrac {\partial z}{\partial \theta }}\\\end{bmatrix}}={\begin{bmatrix}\sin \phi \cos \theta &r\cos \phi \cos \theta &-r\sin \phi \sin \theta \\\sin \phi \sin \theta &r\cos \phi \sin \theta &r\sin \phi \cos \theta \\\cos \phi &-r\sin \phi &0\end{bmatrix}}.$

我們將使用Outer[]函式非常容易地構建這個矩陣，該函式接受除第一個元素之外每個元素的所有可能組合，並將它們作為引數傳遞給第一個元素中的函式。同樣，我們現在不能計算它，因為如果我們這樣做了，Mathematica 將把輸入f_和x_視為“原子”，這意味著它們無法分解，因此不適合作為Outer[].

jacMat[f_, x_] := Outer[D, f, x]

此語句將接受一個方向列表（即 *r*，*φ*，*θ*），以及一個相同長度的函式列表，這些函式是將這些方向對映到目標空間的對映，在本例中是笛卡爾座標系。這個簡單的語句不會檢查你是否給它提供了列表作為輸入，也不會確保它們具有相同的長度。但是，我們知道長度是正確的，所以我們不會為此費心。

現在我們定義方向列表。首先，我們必須移除我們之前對phi和theta的定義（以及r為了安全起見），因為目前，它們將在上述語句中注入一個隨機數，而不是一個符號量。我們還定義函式列表，可以透過將我們現在的符號方向作為引數傳遞給sphericalToCartesian.

sphericalToCartesian[{r_, phi_, theta_}] := {r*Sin[phi]*Cos[theta],  r*Sin[phi]*Sin[theta], r*Cos[phi]};

Remove[r, phi, theta]
dirs = {r, phi, theta}
map = sphericalToCartesian[{r, phi, theta}]

為了找到雅可比矩陣和行列式，我們現在輸入

jacMat[map, dirs]
Simplify[Det[%]]

輸出結果是

{{Cos[theta] Sin[phi], -r Sin[phi] Sin[theta],  r Cos[phi] Cos[theta]}, 
 {Sin[phi] Sin[theta],  r Cos[theta] Sin[phi],  r Cos[phi] Sin[theta]}, 
 {Cos[phi],             0,                     -r Sin[phi]}}

r^2 Sin[phi]

我們看到雅可比矩陣被重新排列了，但在雅可比矩陣中，它是完全等效的，因為它只取決於方向和函式的指定順序。雅可比行列式獨立於經度theta，因此我們在球座標系中的均勻分佈在笛卡爾空間中將是均勻的。但是phi則不會，因為它會隨著雅可比矩陣的變化而變化。r也會這樣做，但由於我們只需要處理一個單位球體，因此我們現在可以將其從等式中刪除。

有一種稱為“逆 CDF 方法”的方法，它允許我們從 0 到 1 之間的一組均勻分佈的數字中生成給定機率密度函式的樣本（就像 Mathematica 的Random[]).

要使用此方法，我們首先必須找到由雅可比行列式給出的 PDF 的 CDF。我們透過除以曲線下的面積來標準化（以給出總機率為 1）。 <>t是一個虛擬變數。

area = Integrate[Sin[t], {t, 0, Pi}];(*Area under PDF*)
cdf[y_] = Integrate[Sin[t]/area, {t, 0, y}](*CDF*)

現在我們找到 CDF 的逆函式

Solve[x == cdf[y], y](*Get distribution transform*)

{{y -> -2 ArcSin[Sqrt[x]]}, {y -> 2 ArcSin[Sqrt[x]]}}

我們取正解，現在我們可以設定我們的theta和phi分佈

 theta := RandomReal[ {0, Pi} ]; (* Uniform on [0,Pi], longitude, only goes halfway round*)
 phi   := 2 ArcSin[Sqrt[Random[]]]; (* Colatitude, from pole to pole *)

現在我們像以前一樣繪製點

dataPoints = Table[{phi, theta}, {6000}];
rectPlot = ListPlot[
  dataPoints,
  AxesLabel -> {"\[Phi]", "\[Theta]"},
  Ticks -> {
    {0, Pi/4, Pi/2, 3 Pi/4, Pi},
    {0, Pi/4, Pi/2, 3 Pi/4, Pi}}
  ]
sphericalToCartesian[{phi_, theta_}] := 
  {Sin[phi]*Cos[theta], 
   Sin[phi]*Sin[theta], 
   Cos[phi]};

mappedPoints = sphericalToCartesian /@ dataPoints;
points = Graphics3D[Point /@ mappedPoints];
Show[points]
Show[points, ViewPoint -> {0, 0, Infinity}]
Show[points, ViewPoint -> {0, Infinity, 0}]

我們減少了極點周圍的分佈密度。

我們現在已經使點在球形區域上均勻分佈。