3D гауссов генератор/преобразование

Простой способ:

Маловероятно, что в такой игре вам понадобятся 3 вариации Гаусса только один раз.

Вам нужна некоторая хранимая переменная, которая может содержать либо триплет гауссовых переменных, либо ничего (Null, Nothing, Empty, что бы это ни было на вашем языке программирования, вы не сказали нам, какое именно).

Изначально в магазине ничего нет (пусто).

Когда попросили тройку:

• если в хранилище содержится триплет, просто верните этот триплет. И пометить магазин как пустой.

• если хранилище пусто, запустите Box-Muller 3 раза. Это дает вам 2 тройки. Ставим вторую тройку в магазин. Верните первую тройку.

Альтернативный способ для математически склонного программиста:

Если кто-то просто попытается адаптировать Box-Muller к 3 измерениям, единственной сложной частью будет получение норма случайного 3D-вектора. Остальное касается двух сферических углов θ (тета) и φ (фи), что несложно.

Оказывается, в 3-х измерениях эта норма включает обратную сторону неполная гамма-функция.

А если у вас Python и Numpy/Scipy, то это функция scipy.special.gammaincinv.

Таким образом, мы можем написать этот код:

import  math
import  numpy.random   as  rd
import  scipy.special  as  sp

# convert 3 uniform [0,1) variates into 3 unit Gaussian variates:
def boxMuller3d(u3):
    u0,u1,u2 = u3    # 3 uniform random numbers in [0,1)

    gamma = u0
    norm2 = 2.0 * sp.gammaincinv(1.5, gamma)  # "regularized" versions
    norm  = math.sqrt(norm2)
    zr    = (2.0 * u1) - 1.0         # sin(theta)
    hr    = math.sqrt(1.0 - zr*zr)   # cos(theta)
    phi   = 2.0 * math.pi * u2
    xr    = hr * math.cos(phi)
    yr    = hr * math.sin(phi)

    g3 = list(map(lambda c: c*norm, [xr, yr, zr]))
    return g3

# generate 3 uniform variates and convert them into 3 unit Gaussian variates:
def gauss3(rng):
    u3 = rng.uniform(0.0, 1.0, 3)
    g3 = boxMuller3d(u3)
    return g3

Чтобы (частично) проверить правильность, у нас может быть эта небольшая основная программа, которая отображает статистические моменты порядка 1-4 результирующего случайного ряда:

randomSeed = 42
rng = rd.default_rng(randomSeed)

count = 3000000  #  (X,Y,Z) triplet count
variates = []
for i  in  range(count):
    g3 = gauss3(rng)
    variates += g3

ln = len(variates)
print("length=%d\n" % ln)

# Checking statistical moments of order 1 to 4:

m1 = sum(variates) / ln
m2 = sum( map(lambda x: x*x,   variates) ) / ln
m3 = sum( map(lambda x: x**3,  variates) ) / ln
m4 = sum( map(lambda x: x**4,  variates) ) / ln

print("m1=%g  m2=%g  m3=%g  m4=%g\n" % (m1,m2,m3,m4))

Вывод тестовой программы:

length=9000000
m1=-0.000455911  m2=1.00025  m3=-0.000563454  m4=3.00184

Таким образом, мы можем видеть, что эти моменты достаточно близки к их математически ожидаемым значениям, соответственно 0,1,0,3.