如何使用 Python 从 3-D 椭球生成点的随机样本？

Question

我试图从 3-D 椭球体中均匀采样大约 1000 个点。有没有某种方法可以对其进行编码，以便我们可以从椭球方程开始获得点？

我想要椭球体表面上的点。

Answer 1

理论

使用这个对 MSE 问题的出色答案如何生成均匀分布在椭球表面上的点？我们可以

在球体上均匀生成一个点，应用映射 f : (x,y,z) -> (x'=ax,y'=by,z'=cz)，然后通过丢弃以一定概率 p(x,y,z) 随机点。

假设椭球体的 3 个轴命名为

0 < a < b < c

我们丢弃生成的点

p(x,y,z) = 1 - mu(x,y,y)/mu_max

概率，即我们将其保留为mu(x,y,y)/mu_max概率，其中

mu(x,y,z) = ((acy)^2 + (abz)^2 + (bcx)^2)^0.5

和

mu_max = bc

执行

import numpy as np
np.random.seed(42)

# Function to generate a random point on a uniform sphere
# (relying on /sf/answers/2378427131/)

def randompoint(ndim=3):
    vec = np.random.randn(ndim,1)
    vec /= np.linalg.norm(vec, axis=0)
    return vec

# Give the length of each axis (example values):

a, b, c = 1, 2, 4

# Function to scale up generated points using the function `f` mentioned above:

f = lambda x,y,z : np.multiply(np.array([a,b,c]),np.array([x,y,z]))

# Keep the point with probability `mu(x,y,z)/mu_max`, ie

def keep(x, y, z, a=a, b=b, c=c):
    mu_xyz = ((a * c * y) ** 2 + (a * b * z) ** 2 + (b * c * x) ** 2) ** 0.5
    return mu_xyz / (b * c) > np.random.uniform(low=0.0, high=1.0)

# Generate points until we have, let's say, 1000 points:

n = 1000
points = []
while len(points) < n:
    [x], [y], [z] = randompoint()
    if keep(x, y, z):
        points.append(f(x, y, z))

支票

检查生成的所有点是否满足椭球条件（即x^2/a^2 + y^2/b^2 + z^2/c^2 = 1）：

for p in points:
    pscaled = np.multiply(p,np.array([1/a,1/b,1/c]))
    assert np.allclose(np.sum(np.dot(pscaled,pscaled)),1)

运行不会出现任何错误。可视化点：

import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(projection="3d")
points = np.array(points)
ax.scatter(points[:, 0], points[:, 1], points[:, 2])
# set aspect ratio for the axes using /sf/answers/4511736281/
ax.set_box_aspect((np.ptp(points[:, 0]), np.ptp(points[:, 1]), np.ptp(points[:, 2])))
plt.show()

这些点看起来分布均匀。

当前接受的答案有问题

在球体上生成一个点，然后仅将其重新投影而不对椭圆进行任何进一步校正，将导致分布扭曲。这本质上与将此帖子 设置为 0 相同。p(x,y,z)想象一个椭球体，其中一个轴比另一个轴大几个数量级。这样，很容易看出，单纯的重投影是行不通的。

Answer 2

这是一个通用函数，用于在具有 a、b 和 c 参数的球体、椭球体或任何三轴椭球体的表面上选取随机点。请注意，直接生成角度不会提供均匀分布，并且会导致沿 z 方向的点数量过多。相反，phi 是作为随机生成的 cos(phi) 的倒数获得的。

    import numpy as np
    def random_point_ellipsoid(a,b,c):
        u = np.random.rand()
        v = np.random.rand()
        theta = u * 2.0 * np.pi
        phi = np.arccos(2.0 * v - 1.0)
        sinTheta = np.sin(theta);
        cosTheta = np.cos(theta);
        sinPhi = np.sin(phi);
        cosPhi = np.cos(phi);
        rx = a * sinPhi * cosTheta;
        ry = b * sinPhi * sinTheta;
        rz = c * cosPhi;
        return rx, ry, rz

该函数取自这篇文章： https: //karthikkaranth.me/blog/generating-random-points-in-a-sphere/