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}]

我们减少了极点周围的分布密度。

我们现在已经使点在球形区域上均匀分布。