从截断的高斯分布中生成 numpy 向量化值
numpy vectorise value generation from a truncated gaussian distribution
我有一个函数，它从 t运行cated 正态分布生成一个值，while 循环确保任何位于 t运行cation 之外的生成值都被丢弃和替换与另一代直到它位于范围内。
def gen_truncated(minimum, maximum, ave, sigma):
# min=0.9, max=1,
x = 0.
while x < minimum or x > maximum:
x = np.random.normal(0,1)*sigma+ave
return x
如何向量化此函数，使 x 现在是许多 x 值的数组，生成方式总是有一个 while 循环确保数组元素只要达到x < minimum和x > maximum的条件就重新生成？是否有将 x 的每个元素与数字进行比较的向量化方法，即 minimum 或 maximum?
编辑：
如果我有更多的约束需要满足怎么办？最终，我希望对通过多个约束生成的 4x4 矩阵的生成进行矢量化，gen_truncated() 中的约束只是众多约束之一。我有一个 gen_sigma()，它首先生成 3 个值 lambda1, lambda2, lambda3，现在 lambda3 再次需要满足几个条件 w.r.t lambda1 和 lambda2 的值否则他们会被重绘。一旦它们正确，所有三个值都会被输入 get_tau() 以生成 3 个值。同样，这些 tau 值需要满足更多约束条件，否则它们将被丢弃并重新生成，直到它们正确为止。最终，它们形成一个名为 sigma_gen 的 4x4 矩阵，通过 gen_channel 左右乘以 create_rotation() 以创建单个 4x4 矩阵 channel.
import numpy as np
from numpy.linalg import norm
def gen_sigma(minimum, maximum, ave, sigma):
lambda1 = gen_truncated(minimum, maximum, ave, sigma)
lambda2 = gen_truncated(minimum, maximum, ave, sigma)
lambda3 = gen_truncated(minimum, maximum, ave, sigma)
while 1+lambda3 < abs(lambda1+lambda2) or 1-lambda3 < abs(lambda2-lambda1):
lambda3 = gen_truncated(minimum, maximum, ave, sigma)
tau = get_tau(lambda1, lambda2, lambda3)
lambdas = [lambda1, lambda2, lambda3]
while (norm(tau)**2 >
1-sum([x**2 for x in [lambda1, lambda2, lambda3]]) +
2*lambda1*lambda2*lambda3) or (z_eta(tau, lambdas) < 0):
tau = get_tau(lambda1, lambda2, lambda3)
sigma_gen = np.array([[ 1, 0, 0, 0],
[tau[0], lambda1, 0, 0],
[tau[1], 0, lambda2, 0],
[tau[2], 0, 0, lambda3]])
return sigma_gen
def get_tau(einval1, einval2, einval3):
max_tau1 = 1 - abs(einval1)
max_tau2 = 1 - abs(einval2)
max_tau3 = 1 - abs(einval3)
tau1 = max_tau1*(2*np.random.uniform(0,1)-1)
tau2 = max_tau2*(2*np.random.uniform(0,1)-1)
tau3 = max_tau3*(2*np.random.uniform(0,1)-1)
return [tau1, tau2, tau3]
def z_eta(t: np.ndarray, l: np.ndarray):
condition = (norm(t)**4 - 2*norm(t)**2 -
2*sum([(l[i]**2)*(2*(t[i]**2-norm(t)**2)) for i in range(3)])+
q(l))
return condition
def q(e: np.ndarray):
# e are the eigenvalues
return (1+e[0]+e[1]+e[2])*(1+e[0]-e[1]-e[2])*(1-e[0]+e[1]-e[2])*(1-e[0]-e[1]+e[2])
def create_rotation(angles: np.ndarray) -> np.ndarray:
"random rotation in PL form"
# input np.random.normal(0,1,3)*0.06
rotation = np.eye(4, dtype=complex)
left = np.array([[ np.cos(angles[0]), np.sin(angles[0]), 0],
[-np.sin(angles[0]), np.cos(angles[0]), 0],
[ 0, 0, 1]])
mid = np.array([[1, 0, 0],
[0, np.cos(angles[1]), np.sin(angles[1])],
[0, -np.sin(angles[1]), np.cos(angles[1])]])
right = np.array([[ np.cos(angles[2]), np.sin(angles[2]), 0],
[-np.sin(angles[2]), np.cos(angles[2]), 0],
[ 0, 0, 1]])
rotation[1:4,1:4] = left@mid@right
return rotation
def gen_channel(r1, r2, ave, sigma):
rand1 = np.random.normal(0,1,3)
rand2 = np.random.normal(0,1,3)
channel = create_rotation(rand1*r1)@gen_sigma(0.9, 1, ave, sigma)@\
create_rotation(rand2*r2)
return channel
频道示例运行
gen_channel(0.05, 0.05, 0.98, 0.15)
例如
Out[140]:
array([[ 1. +0.j, 0. +0.j, 0. +0.j,
0. +0.j],
[-0.05828008+0.j, 0.91805971+0.j, 0.14291751+0.j,
-0.00946994+0.j],
[-0.00509449+0.j, -0.14170308+0.j, 0.90034613+0.j,
-0.11548884+0.j],
[ 0.0467522 +0.j, -0.00851749+0.j, 0.11450963+0.j,
0.90259637+0.j]])
现在，如果我想创建这些 4x4 矩阵中的 100 个，我将不得不使用列表理解，即
np.array([gen_channel(0.05, 0.05, 0.98, 0.15) for i in range(100)])
这将 运行 通过所有约束比较并一一创建 4x4 矩阵。现在我最初的问题是因为我想对它们进行矢量化，所以与其一次比较一个值，不如使用 numpy broadcast 生成一组值并检查约束，这样我就有了 [= 的矢量化版本31=] 生成 100 个这样的 4x4 矩阵而不需要列表理解。列表理解方式包含重复使用生成单个随机数，这导致其运行速度出现瓶颈。我想做的只是生成随机数数组，进行这些检查，然后生成 4x4 通道数组以减少瓶颈。
您可以从原始分布中抽取大量样本，然后确定哪些条目位于正确的范围内，然后从中抽取：
# parameters
ave, sigma = 0,1
minimum, maximum = 0.9, 1
# draw sample and specify which entries are ok
a = np.random.normal(ave, sigma, 100000)
index = (a > minimum) & (a < maximum)
# draw from subset
np.random.choice(a[index], 1000, replace=False)
使用 timeit
关于上面的代码：
%%timeit -r 10 -n 10
2.51 ms ± 87.5 µs per loop (mean ± std. dev. of 10 runs, 10 loops each)
在循环中的原件上：
%%timeit -r 10 -n 10
for i in range(1000):
gen_truncated(0.9,1, 0, 1)
88.5 ms ± 1.24 ms per loop (mean ± std. dev. of 10 runs, 10 loops each)
我有一个函数，它从 t运行cated 正态分布生成一个值，while 循环确保任何位于 t运行cation 之外的生成值都被丢弃和替换与另一代直到它位于范围内。
def gen_truncated(minimum, maximum, ave, sigma):
# min=0.9, max=1,
x = 0.
while x < minimum or x > maximum:
x = np.random.normal(0,1)*sigma+ave
return x
如何向量化此函数，使 x 现在是许多 x 值的数组，生成方式总是有一个 while 循环确保数组元素只要达到x < minimum和x > maximum的条件就重新生成？是否有将 x 的每个元素与数字进行比较的向量化方法，即 minimum 或 maximum?
编辑：
如果我有更多的约束需要满足怎么办？最终，我希望对通过多个约束生成的 4x4 矩阵的生成进行矢量化，gen_truncated() 中的约束只是众多约束之一。我有一个 gen_sigma()，它首先生成 3 个值 lambda1, lambda2, lambda3，现在 lambda3 再次需要满足几个条件 w.r.t lambda1 和 lambda2 的值否则他们会被重绘。一旦它们正确，所有三个值都会被输入 get_tau() 以生成 3 个值。同样，这些 tau 值需要满足更多约束条件，否则它们将被丢弃并重新生成，直到它们正确为止。最终，它们形成一个名为 sigma_gen 的 4x4 矩阵，通过 gen_channel 左右乘以 create_rotation() 以创建单个 4x4 矩阵 channel.
import numpy as np
from numpy.linalg import norm
def gen_sigma(minimum, maximum, ave, sigma):
lambda1 = gen_truncated(minimum, maximum, ave, sigma)
lambda2 = gen_truncated(minimum, maximum, ave, sigma)
lambda3 = gen_truncated(minimum, maximum, ave, sigma)
while 1+lambda3 < abs(lambda1+lambda2) or 1-lambda3 < abs(lambda2-lambda1):
lambda3 = gen_truncated(minimum, maximum, ave, sigma)
tau = get_tau(lambda1, lambda2, lambda3)
lambdas = [lambda1, lambda2, lambda3]
while (norm(tau)**2 >
1-sum([x**2 for x in [lambda1, lambda2, lambda3]]) +
2*lambda1*lambda2*lambda3) or (z_eta(tau, lambdas) < 0):
tau = get_tau(lambda1, lambda2, lambda3)
sigma_gen = np.array([[ 1, 0, 0, 0],
[tau[0], lambda1, 0, 0],
[tau[1], 0, lambda2, 0],
[tau[2], 0, 0, lambda3]])
return sigma_gen
def get_tau(einval1, einval2, einval3):
max_tau1 = 1 - abs(einval1)
max_tau2 = 1 - abs(einval2)
max_tau3 = 1 - abs(einval3)
tau1 = max_tau1*(2*np.random.uniform(0,1)-1)
tau2 = max_tau2*(2*np.random.uniform(0,1)-1)
tau3 = max_tau3*(2*np.random.uniform(0,1)-1)
return [tau1, tau2, tau3]
def z_eta(t: np.ndarray, l: np.ndarray):
condition = (norm(t)**4 - 2*norm(t)**2 -
2*sum([(l[i]**2)*(2*(t[i]**2-norm(t)**2)) for i in range(3)])+
q(l))
return condition
def q(e: np.ndarray):
# e are the eigenvalues
return (1+e[0]+e[1]+e[2])*(1+e[0]-e[1]-e[2])*(1-e[0]+e[1]-e[2])*(1-e[0]-e[1]+e[2])
def create_rotation(angles: np.ndarray) -> np.ndarray:
"random rotation in PL form"
# input np.random.normal(0,1,3)*0.06
rotation = np.eye(4, dtype=complex)
left = np.array([[ np.cos(angles[0]), np.sin(angles[0]), 0],
[-np.sin(angles[0]), np.cos(angles[0]), 0],
[ 0, 0, 1]])
mid = np.array([[1, 0, 0],
[0, np.cos(angles[1]), np.sin(angles[1])],
[0, -np.sin(angles[1]), np.cos(angles[1])]])
right = np.array([[ np.cos(angles[2]), np.sin(angles[2]), 0],
[-np.sin(angles[2]), np.cos(angles[2]), 0],
[ 0, 0, 1]])
rotation[1:4,1:4] = left@mid@right
return rotation
def gen_channel(r1, r2, ave, sigma):
rand1 = np.random.normal(0,1,3)
rand2 = np.random.normal(0,1,3)
channel = create_rotation(rand1*r1)@gen_sigma(0.9, 1, ave, sigma)@\
create_rotation(rand2*r2)
return channel
频道示例运行
gen_channel(0.05, 0.05, 0.98, 0.15)
例如
Out[140]:
array([[ 1. +0.j, 0. +0.j, 0. +0.j,
0. +0.j],
[-0.05828008+0.j, 0.91805971+0.j, 0.14291751+0.j,
-0.00946994+0.j],
[-0.00509449+0.j, -0.14170308+0.j, 0.90034613+0.j,
-0.11548884+0.j],
[ 0.0467522 +0.j, -0.00851749+0.j, 0.11450963+0.j,
0.90259637+0.j]])
现在，如果我想创建这些 4x4 矩阵中的 100 个，我将不得不使用列表理解，即
np.array([gen_channel(0.05, 0.05, 0.98, 0.15) for i in range(100)])
这将 运行 通过所有约束比较并一一创建 4x4 矩阵。现在我最初的问题是因为我想对它们进行矢量化，所以与其一次比较一个值，不如使用 numpy broadcast 生成一组值并检查约束，这样我就有了 [= 的矢量化版本31=] 生成 100 个这样的 4x4 矩阵而不需要列表理解。列表理解方式包含重复使用生成单个随机数，这导致其运行速度出现瓶颈。我想做的只是生成随机数数组，进行这些检查，然后生成 4x4 通道数组以减少瓶颈。
您可以从原始分布中抽取大量样本，然后确定哪些条目位于正确的范围内，然后从中抽取：
# parameters
ave, sigma = 0,1
minimum, maximum = 0.9, 1
# draw sample and specify which entries are ok
a = np.random.normal(ave, sigma, 100000)
index = (a > minimum) & (a < maximum)
# draw from subset
np.random.choice(a[index], 1000, replace=False)
使用 timeit
关于上面的代码：
%%timeit -r 10 -n 10
2.51 ms ± 87.5 µs per loop (mean ± std. dev. of 10 runs, 10 loops each)
在循环中的原件上：
%%timeit -r 10 -n 10
for i in range(1000):
gen_truncated(0.9,1, 0, 1)
88.5 ms ± 1.24 ms per loop (mean ± std. dev. of 10 runs, 10 loops each)