你是这样做的，我用10代替了10000，因为这需要一段时间。我最初的猜测是0，我把它设置为上一次迭代的下一个猜测，因为它应该非常接近解决方案。如果你愿意，你可以进一步限制它，所以它严格高于它。
作为一个侧面的评论，这种复杂分布的抽样是不可行的，因为计算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 = 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是如何工作的。