离散正态分布以获得随机变量的概率

Question

假设我从一个均值为零的正态分布中随机抽取，标准差由一个向量表示，例如，维度为 3 的向量

scale_rng=np.array([1,2,3])
eps=np.random.normal(0,scale_rng)

我需要根据我绘制上述eps 的一些模拟计算加权平均值。这个平均值的权重是“eps 的概率”（因此我将有一个具有 3 个权重的向量）。对于加权平均，我只是指算术和，其中每个分量乘以权重，即 0 到 1 之间的数字，并且所有权重的总和应为 1。 这种加权平均值应计算如下：我有一个变量的时间序列观察，x。我计算了x 的扩展滚动标准偏差（假设这是scale 中的值）。然后，我从正态分布中提取一个随机变量eps，如上所述，对于x 中的每个时间观察，我将其添加到其中，比如获得y=x+eps。最后，我需要计算y 的加权平均值，其中y 的每个值由“从均值为零且标准差等于scale 的正态分布中提取eps 的每个值的概率加权。

现在，我知道我不能认为这是 pdf 上与随机绘制的值相对应的点，因为正常的随机变量是连续的，因此某个点的 pdf 为零。因此，我发现的唯一解决方案是用一定数量的 bin 离散正态分布，然后找到用上面的代码提取的值实际被绘制的概率。我怎么能在 Python 中做到这一点？

编辑：我找到的解决方案是使用

norm.cdf(eps_it+0.5, loc=0, scale=scale_rng)-norm.cdf(eps_it-0.5, loc=0, scale=scale_rng)

这并不是真正基于离散化，但至少在我看来“概率方面”是可行的。

Answer 1

这是一个让所有内容保持连续的示例。

import numpy as np
from scipy import stats

# some function we want a monte carlo estimate of
def fn(eps):
  return np.sum(np.abs(eps), axis=1)

# define distribution of eps
sd = np.array([1,2,3])
d_eps = stats.norm(0, sd)

# draw uniform samples so we don't double apply the normal density
eps = np.random.uniform(-6*sd, 6*sd, size=(10000, 3))

# calculate weights (working with log-likelihood is better for numerical stability)
w = np.prod(d_eps.pdf(eps), axis=1)
# normalise so weights sum to 1
w /= np.sum(w)

# get estimate
np.sum(fn(eps) * w)

如果我运行几次，它会给我 4.71、4.74、4.70 4.78。当eps 直接从法线中提取时，我们可以通过使用平均值来验证这是正确的：

np.mean(fn(d_eps.rvs(size=(10000, 3))))

这给了我基本相同的值，但预期的方差更低。例如4.79、4.76、4.77、4.82、4.80。