将Gaussian分布转换为总粒子数为N的整数型1D直方图的方法求助

实现方案 #

核心逻辑 #

你需要解决两个核心问题：将高斯分布的概率权重归一化到总粒子数N、保证每个bin的粒子数为整数且总和严格等于N，可以按以下步骤实现：

• 计算每个网格bin对应的高斯概率权重，你原来的高斯函数输出的是概率密度，乘以bin宽度（这里为1，可省略）得到单个bin的采样概率
• 将所有bin的概率归一化后乘以N，得到每个bin的预期粒子数（浮点型）
• 先取每个预期值的整数部分作为基础粒子数，统计未分配的剩余粒子数
• 将剩余粒子按每个bin的余数从大到小排序，依次给每个bin分配1个粒子，直到剩余粒子全部分配完成，该方案可以保证最终分布和理论高斯分布的偏差最小

完整可运行代码 #

import numpy as np

def gaussian(x, mu, sig):
    return 1./(np.sqrt(2.*np.pi)*sig)*np.exp(-np.power((x - mu)/sig, 2.)/2)

# 配置参数
N = 1000  # 总粒子数，可自行修改
mean = 0
sigma = 1
x_min = -10
x_max = 10
bin_width = 1

# 生成网格点
x_values = np.arange(x_min, x_max, bin_width)
# 计算每个bin的概率密度
pdf_values = gaussian(x_values, mean, sigma)
# 归一化得到权重，乘以N得到浮点预期粒子数
expected_float = pdf_values / pdf_values.sum() * N

# 取整数部分作为基础粒子数
counts = np.floor(expected_float).astype(int)
# 计算剩余需要分配的粒子数
remainder = expected_float - counts
remaining_particles = int(N - counts.sum())

# 按余数从大到小排序的索引，取前remaining_particles个各加1
top_remainder_indices = remainder.argsort()[-remaining_particles:]
counts[top_remainder_indices] += 1

# 输出验证
print(f"总粒子数：{counts.sum()}，等于设定的N={N}")
print(f"各位置粒子数：\n{list(zip(x_values, counts))}")

结果验证 #

• 最终输出的counts数组就是每个网格位置对应的粒子数，总和严格等于N
• 粒子分布完全符合设定均值和标准差的高斯分布，不存在你之前遇到的平直分布问题
• 如果需要生成每个粒子的具体位置坐标而不是直方图计数，可以直接用np.repeat(x_values, counts)得到长度为N的粒子位置数组

内容的提问来源于stack exchange，提问作者ValientProcess