Python：使用逆 cdf 技术进行采样答案

【问题标题】：Python: Sampling using inverse cdf techniquePython：使用逆 cdf 技术进行采样
【发布时间】：2015-12-16 17:01:50
【问题描述】：

我有一个复杂的（非标准）分布函数，我想对其进行采样以使用逆 cdf 技术生成模拟数据点。 为了这个例子，我将考虑一个高斯分布

var=100
def f(x,a):
     def g(y):
       return (1/np.sqrt(2*np.pi*var))*np.exp(-y**2/(2*var))
   b,err=integrate.quad(g,-np.inf,x) 
   return b-a

我想在a=[0,1]、a=np.linspace(0,1,10000,endpoint=False) 之间生成值，并使用scipy.optimize.fsolve 为每个a 求解x。我有两个问题：

如何将fsolve 用于值数组a？
fsolve 进行初始猜测x0，如何选择一个好的猜测值？

谢谢

【问题讨论】：

这里a 的目的是什么，您是否要为a 中的每个值反转cdf？我认为fsolve 不会以这种方式处理数组，您可能需要调用它 10000 次。
关于 2 我不太确定，但您可以根据之前的 a 绑定解决方案，因为 cdf 是非递减的。
@simonzack a 是 icdf 的值，我在定义中减去它让 fsolve 求解方程。不使用 fsolve 有没有其他数值求解积分上限的方法？
我不认为这个问题是解决它的地方，但为了提高效率，我建议使用ODE integrator 来集成 PDF once，然后将将结果值转换为插值对象以创建逆 CDF。之后，评估逆 CDF 只是一个函数调用，它是矢量化的（它适用于数组...我认为），因此您可以像 icdf(np.random.uniform(size=10000)) 一样轻松获取样本。

标签： python scipy

【解决方案1】：

你是这样做的，我用 10 替换了 10000，因为这需要一段时间。我最初的猜测只是 0，我将它设置为上一次迭代以进行下一次猜测，因为它应该非常接近解决方案。如果你愿意，你可以进一步绑定它，所以它严格高于它。

作为旁注，这种复杂分布的采样实际上并不可行，因为计算 cdf 可能相当困难。还有其他采样技术可以解决这些问题，例如 Gibbs 采样、Metropolis Hastings 等。

var = 100

def f(x, a):
    def g(y):
        return (1/np.sqrt(2*np.pi*var))*np.exp(-y**2/(2*var))

    b, err = sp.integrate.quad(g, -np.inf, x) 
    return b - a


a = np.linspace(0, 1, 10, endpoint=False)[1:]
x0 = 0
for a_ in a:
    xi = sp.optimize.fsolve(f, x0 + 0.01, args=(a_,))[0]
    print(xi)
    x0 = xi

[编辑] 它似乎卡在 0 附近，添加一个小数字可以解决它，我不知道为什么，因为我不知道 fsolve 是如何工作的。

【讨论】：

我认为this question 谈到了卡住的原因。
@DavidZ 这确实有些道理。
谢谢，您的解决方案有效，我选择初始猜测为 1 而不是 0。至于您对使用 Metropolis-Hasting 和其他 MCMC 技术的评论，是的，我确实尝试使用 emcee 来采样分布，但是当我使用模拟值来最大化似然性时，我的最佳拟合值已经偏离了！也许我会打开一个新帖子，重点介绍我所做的事情，看看是否有办法提高准确性。