生成随机数a、b、c,使a^2+b^2+c^2=1

根据 this 回答统计信息。stackexchange。com,您应该使用正态分布值来获得球体上均匀分布的值。这意味着,你可以:

import numpy as np
abc = np.random.normal(size=3)
a,b,c = abc/np.sqrt(sum(abc**2))

import numpy as np
import random
import math

def MSeifert():
    a = 1
    b = 1
    while a**2 + b**2 > 1:  # discard any a and b whose sum of squares already exceeds 1
        a = random.random()
        b = random.random()
    c = math.sqrt(1 - a**2 - b**2)  # fixed c
    return a, b, c

def VBB():
    x = np.random.uniform(0,1,3) # random numbers in [0, 1)
    x /= np.sqrt(x[0] ** 2 + x[1] ** 2 + x[2] ** 2)
    return x[0], x[1], x[2]

def user3684792():
    theta = random.uniform(0, 0.5*np.pi)
    phi = random.uniform(0, 0.5*np.pi)
    return np.sin(theta)* np.cos(phi), np.sin(theta)*np.sin(phi), np.cos(theta)

def JohanL():
    abc = np.random.normal(size=3)
    a,b,c = abc/np.sqrt(sum(abc**2))
    return a, b, c

def SeverinPappadeux():
    cos_th = 2.0*random.uniform(0, 1.0) - 1.0
    sin_th = math.sqrt(1.0 - cos_th*cos_th)
    phi = random.uniform(0, 2.0*math.pi)
    return sin_th * math.cos(phi), sin_th * math.sin(phi), cos_th

并绘制分布图:

%matplotlib notebook

import matplotlib.pyplot as plt

f, axes = plt.subplots(3, 4)

for func_idx, func in enumerate([MSeifert, JohanL, user3684792, VBB]):
    axes[0, func_idx].set_title(str(func.__name__))
    res = [func() for _ in range(50000)]
    for idx in range(3):
        axes[idx, func_idx].hist([i[idx] for i in res], bins='auto')

axes[0, 0].set_ylabel('a')
axes[1, 0].set_ylabel('b')
axes[2, 0].set_ylabel('c')

plt.tight_layout()

说明:这些行显示了 a , b 和 c

在该范围内给出均匀随机分布的唯一方法 (-1, 1) 是约翰和塞维林·帕帕德的方法。所有其他方法都有一些特征,如尖峰或范围内的功能行为 [0, 1) . 注意,这两种解决方案目前给出的值介于-1和1之间,而所有其他方法给出的值介于0和1之间。

我认为这实际上是一个很酷的问题,一个很好的方法是只使用球面极坐标并随机生成角度。

import random
import numpy as np
def random_pt():
    theta = random.uniform(0, 0.5*np.pi)
    phi = random.uniform(0, 0.5*np.pi)
    return np.sin(theta)* np.cos(phi), np.sin(theta)*np.sin(phi), np.cos(theta)

你可以这样做:

import random
import math

def three_random_numbers_adding_to_one():
    a = 1
    b = 1
    while a**2 + b**2 > 1:  # discard any a and b whose sum of squares already exceeds 1
        a = random.random()
        b = random.random()
    c = math.sqrt(1 - a**2 - b**2)  # fixed c
    return a, b, c

a, b, c = three_random_numbers_adding_to_one()
print(a**2 + b**2 + c**2)

然而 float 1 ,大约。

您可能需要检查使用此函数生成的数字是否“足够随机”。这可能是因为这种设置偏向了“随机性”。

“正确”的答案取决于您是在空间中寻找均匀的随机分布,还是在球体的表面上,或者其他什么。如果要在球体表面上寻找点,仍然需要担心 cos(theta) 导致点在球体极点附近“聚在一起”的因素。由于从你的问题中不清楚确切的性质,这里有一个“完全随机”的分布应该有效:

x = np.random.uniform(0,1,3) # random numbers in [0, 1)
x /= np.sqrt(x[0] ** 2 + x[1] ** 2 + x[2] ** 2)

这里的另一个优点是,由于我们使用numpy阵列,您也可以通过使用 x = np.random.uniform(0, 1, (3, n)) 对于任何 n

是时候添加另一个解决方案了,呵呵。。。

这一次,它在单位球体点拾取-检查上是真正一致的 http://mathworld.wolfram.com/SpherePointPicking.html 有关详细信息

import math
import random

def random_pt():
    cos_th = 2.0*random.uniform(0, 1.0) - 1.0
    sin_th = math.sqrt(1.0 - cos_th*cos_th)
    phi = random.uniform(0, 2.0*math.pi)
    return sin_th * math.cos(phi), sin_th * math.sin(phi), cos_th

for k in range(0, 100):
    a, b, c = random_pt()

    print("{0} {1} {2} {3}".format(a, b, c, a*a + b*b + c*c))