生成长度分布为 exp(-x/lambda) 的各向同性步长
Generating isotropic steps with lengths distributed as exp(-x/lambda)
我正在编写有关辐射传输的代码，并已生成各向同性矢量，如我的代码所示。我现在不确定如何生成长度分布为 exp(-x/lambda) 的各向同性步长。有什么建议吗？
def isotropic_unit_vectors():
nparticles = 1000
lambda_a = 45
lambda_s = 0.3
sigma_a = 1/lambda_a
sigma_s = 1/lambda_s
sigma_T = sigma_a + sigma_s
lambda_T = 1/sigma_T
u = np.random.uniform(low = 0, high =1, size = 1000)
step = -lambda_T*np.log(u)
theta = 2*np.pi*np.random.rand(nparticles)
phi = np.arccos(2*np.random.rand(nparticles)-1)
dx = step*np.cos(theta)*np.cos(phi)
dy = step*np.cos(theta)*np.sin(phi)
dz = step*np.sin(theta)
rand = np.column_stack((dx,dy,dz))
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(dx, dy, dz, c='r', marker='.')
ax.set_xlabel('X Label')
ax.set_ylabel('Y Label')
ax.set_zlabel('Z Label')
plt.show()
return rand
编辑：我希望多个粒子从 0,0,0 开始，然后随机移动。
我假设你想用 CDF 生成随机数； F(x) = exp(-x/lambda)，对于 x 的有效支持，使得 0 <= F(x) <= 1.
事实：
任意随机数的 CDF 均匀分布为 uniform(0,1)（CDF 逆方法）。也就是说，生成一个统一的 (0,1) 随机数，并将其设置为等于以 x 表示的 CDF。然后求解 x（参见步骤 2）。
方法
您可以遵循以下流程：
- 在(0, 1)之间均匀生成一个随机数，我们称是u1。
- 设置 u1 = exp(-x/lambda) => x = -lambda*log(u1)
- 将 x 保存在列表或其他内容中。
- 返回步骤1，直到生成所需数量的随机数。
注：
如果 exp(-x/lambda) 是 PDF 而不是 CDF，您可以通过集成 PDF 找到 CDF 并执行上述相同的步骤，但在步骤 2 中有所更改，具体取决于 CDF 是什么。对于第 1 步，您可以使用 random.random() 在 Py 中生成一个 unif(0,1)。
希望这对您有所帮助。
好的，你的点点滴滴都是正确的，但它们放在一起并不完全正确。它可能看起来像，Python 3.7，Anaconda WIndows 10 x64
import numpy as np
def isotropic_unit_vectors():
"""Generates troika of directional cosines isotropically on the sphere"""
cs_theta = 2.0 * np.random.random() - 1.0
sn_theta = np.sqrt((1.0 - cs_theta)*(1.0 + cs_theta))
phi = 2.0 * np.pi * np.random.random()
wx = sn_theta * np.cos(phi)
wy = sn_theta * np.sin(phi)
wz = cs_theta
return (wx, wy, wz)
def random_step(λ):
u = 1.0 - np.random.random() # NB! avoids log(0)
return -λ * np.log(u)
nparticles = 100
λ_a = 45 # absorbtion
λ_s = 0.3 # scattering
# cross-sections
σ_a = 1/λ_a
σ_s = 1/λ_s
σ_t = σ_a + σ_s
λ_t = 1/σ_t
trajectories = {} # dictionary for all trajectories
for k in range(0, nparticles):
trajectory = []
# initialize particle position
x = 0.0
y = 0.0
z = 0.0
print(f"Partile number {k+1}")
print(f" Position {x}, {y}, {z}")
trajectory.append((x, y, z))
while True:
(wx, wy, wz) = isotropic_unit_vectors() # sample direction
step = random_step(λ_t)
# make step and find new position
x += wx * step
y += wy * step
z += wz * step
print(f" Position {x}, {y}, {z}")
trajectory.append((x, y, z))
# simulate absorbtion
if np.random.random() < σ_a/σ_t:
trajectories[k] = trajectory
break # absorbtion, end of trajectory, next particle would be tracked
else:
pass # scattering, continue with this trajectory
请注意，代码未向量化，可以使用 NumPy 向量做得更好更快。
更新
添加了代码以将所有轨迹保存在由 # 索引的字典中。轨迹本身是元组列表，每个元组包含行走的3个坐标
更新二
如果你在得到trajectories字典后添加代码，它会绘制一个带有顶点的轨迹。顶点颜色由距 (0,0,0) 的距离定义。
from mpl_toolkits import mplot3d
import matplotlib.pyplot as plt
fig = plt.figure()
ax = plt.axes(projection='3d')
# Data for a three-dimensional line
k = 7 # print particle #8
trj = trajectories[k]
xline = [x for x, y, z in trj]
yline = [y for x, y, z in trj]
zline = [z for x, y, z in trj]
# plot trajectory
ax.plot3D(xline, yline, zline, 'gray')
# plot vertices
xpts = [x for x, y, z in trj]
ypts = [y for x, y, z in trj]
zpts = [z for x, y, z in trj]
rpts = [np.sqrt(x*x+y*y+z*z) for x, y, z in trj]
ax.scatter3D(xpts, ypts, zpts, c=rpts, cmap='hsv')
plt.show()
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
def isotropic_step_generator(l, size):
# l is the attenuation length (it's the l in the distribution exp(-x/l))
theta = theta = np.arccos(1-2*np.random.uniform(0, 1, size))
phi = np.random.uniform(0, 2*np.pi, size)
radius = -l * np.log(np.random.uniform(size=size))
x = radius * np.sin(theta) * np.cos(phi)
y = radius * np.sin(theta) * np.sin(phi)
z = radius * np.cos(theta)
return x, y, z
empty_array = np.empty((1000,3))
i=0
while i<1000:
empty_array[i] = isotropic_step_generator(45, 1) # attenuation length for water is 45 cm
i+=1
x_array = empty_array[:, 0]
y_array = empty_array[:, 1]
z_array = empty_array[:, 2]
fig = plt.figure()
ax = Axes3D(fig)
ax.scatter(x_array, y_array, z_array, s=7)
plt.show()
以上代码生成长度分布为 exp(-x/lambda) 的各向同性步长，并将它们绘制为水。请注意，您必须对 theta 使用 arccos 以防止点聚集到极点。
我正在编写有关辐射传输的代码，并已生成各向同性矢量，如我的代码所示。我现在不确定如何生成长度分布为 exp(-x/lambda) 的各向同性步长。有什么建议吗？
def isotropic_unit_vectors():
nparticles = 1000
lambda_a = 45
lambda_s = 0.3
sigma_a = 1/lambda_a
sigma_s = 1/lambda_s
sigma_T = sigma_a + sigma_s
lambda_T = 1/sigma_T
u = np.random.uniform(low = 0, high =1, size = 1000)
step = -lambda_T*np.log(u)
theta = 2*np.pi*np.random.rand(nparticles)
phi = np.arccos(2*np.random.rand(nparticles)-1)
dx = step*np.cos(theta)*np.cos(phi)
dy = step*np.cos(theta)*np.sin(phi)
dz = step*np.sin(theta)
rand = np.column_stack((dx,dy,dz))
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(dx, dy, dz, c='r', marker='.')
ax.set_xlabel('X Label')
ax.set_ylabel('Y Label')
ax.set_zlabel('Z Label')
plt.show()
return rand
编辑：我希望多个粒子从 0,0,0 开始，然后随机移动。
我假设你想用 CDF 生成随机数； F(x) = exp(-x/lambda)，对于 x 的有效支持，使得 0 <= F(x) <= 1.
事实： 任意随机数的 CDF 均匀分布为 uniform(0,1)（CDF 逆方法）。也就是说，生成一个统一的 (0,1) 随机数，并将其设置为等于以 x 表示的 CDF。然后求解 x（参见步骤 2）。
方法 您可以遵循以下流程：
- 在(0, 1)之间均匀生成一个随机数，我们称是u1。
- 设置 u1 = exp(-x/lambda) => x = -lambda*log(u1)
- 将 x 保存在列表或其他内容中。
- 返回步骤1，直到生成所需数量的随机数。
注：
如果 exp(-x/lambda) 是 PDF 而不是 CDF，您可以通过集成 PDF 找到 CDF 并执行上述相同的步骤，但在步骤 2 中有所更改，具体取决于 CDF 是什么。对于第 1 步，您可以使用 random.random() 在 Py 中生成一个 unif(0,1)。
希望这对您有所帮助。
好的，你的点点滴滴都是正确的，但它们放在一起并不完全正确。它可能看起来像，Python 3.7，Anaconda WIndows 10 x64
import numpy as np
def isotropic_unit_vectors():
"""Generates troika of directional cosines isotropically on the sphere"""
cs_theta = 2.0 * np.random.random() - 1.0
sn_theta = np.sqrt((1.0 - cs_theta)*(1.0 + cs_theta))
phi = 2.0 * np.pi * np.random.random()
wx = sn_theta * np.cos(phi)
wy = sn_theta * np.sin(phi)
wz = cs_theta
return (wx, wy, wz)
def random_step(λ):
u = 1.0 - np.random.random() # NB! avoids log(0)
return -λ * np.log(u)
nparticles = 100
λ_a = 45 # absorbtion
λ_s = 0.3 # scattering
# cross-sections
σ_a = 1/λ_a
σ_s = 1/λ_s
σ_t = σ_a + σ_s
λ_t = 1/σ_t
trajectories = {} # dictionary for all trajectories
for k in range(0, nparticles):
trajectory = []
# initialize particle position
x = 0.0
y = 0.0
z = 0.0
print(f"Partile number {k+1}")
print(f" Position {x}, {y}, {z}")
trajectory.append((x, y, z))
while True:
(wx, wy, wz) = isotropic_unit_vectors() # sample direction
step = random_step(λ_t)
# make step and find new position
x += wx * step
y += wy * step
z += wz * step
print(f" Position {x}, {y}, {z}")
trajectory.append((x, y, z))
# simulate absorbtion
if np.random.random() < σ_a/σ_t:
trajectories[k] = trajectory
break # absorbtion, end of trajectory, next particle would be tracked
else:
pass # scattering, continue with this trajectory
请注意，代码未向量化，可以使用 NumPy 向量做得更好更快。
更新
添加了代码以将所有轨迹保存在由 # 索引的字典中。轨迹本身是元组列表，每个元组包含行走的3个坐标
更新二
如果你在得到trajectories字典后添加代码，它会绘制一个带有顶点的轨迹。顶点颜色由距 (0,0,0) 的距离定义。
from mpl_toolkits import mplot3d
import matplotlib.pyplot as plt
fig = plt.figure()
ax = plt.axes(projection='3d')
# Data for a three-dimensional line
k = 7 # print particle #8
trj = trajectories[k]
xline = [x for x, y, z in trj]
yline = [y for x, y, z in trj]
zline = [z for x, y, z in trj]
# plot trajectory
ax.plot3D(xline, yline, zline, 'gray')
# plot vertices
xpts = [x for x, y, z in trj]
ypts = [y for x, y, z in trj]
zpts = [z for x, y, z in trj]
rpts = [np.sqrt(x*x+y*y+z*z) for x, y, z in trj]
ax.scatter3D(xpts, ypts, zpts, c=rpts, cmap='hsv')
plt.show()
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
def isotropic_step_generator(l, size):
# l is the attenuation length (it's the l in the distribution exp(-x/l))
theta = theta = np.arccos(1-2*np.random.uniform(0, 1, size))
phi = np.random.uniform(0, 2*np.pi, size)
radius = -l * np.log(np.random.uniform(size=size))
x = radius * np.sin(theta) * np.cos(phi)
y = radius * np.sin(theta) * np.sin(phi)
z = radius * np.cos(theta)
return x, y, z
empty_array = np.empty((1000,3))
i=0
while i<1000:
empty_array[i] = isotropic_step_generator(45, 1) # attenuation length for water is 45 cm
i+=1
x_array = empty_array[:, 0]
y_array = empty_array[:, 1]
z_array = empty_array[:, 2]
fig = plt.figure()
ax = Axes3D(fig)
ax.scatter(x_array, y_array, z_array, s=7)
plt.show()
以上代码生成长度分布为 exp(-x/lambda) 的各向同性步长，并将它们绘制为水。请注意，您必须对 theta 使用 arccos 以防止点聚集到极点。