def expect_mc_bounds(dist, func=lambda x: 1, size=50000, lower=None, upper=None,
                     conditional=False, overfact=1.2):
    #call rvs once to find length of random vector
    rvsdim = dist.rvs(size=1).shape[-1]
    if lower is None:
        lower = -np.inf * np.ones(rvsdim)
    else:
        lower = np.asarray(lower)
    if upper is None:
        upper = np.inf * np.ones(rvsdim)
    else:
        upper = np.asarray(upper)

    def fun(x):
        return func(x) # * dist.pdf(x)

    rvsli = []
    used = 0 #remain = size  #inplace changes size
    total = 0
    while True:
        remain = size - used  #just a temp variable
        rvs = dist.rvs(size=int(remain * overfact))
        total += int(size * overfact)

        rvsok = rvs[((rvs >= lower) & (rvs <= upper)).all(-1)]
        #if rvsok.ndim == 1: #possible shape problems if only 1 random vector
        rvsok = np.atleast_2d(rvsok)
        used += rvsok.shape[0]

        rvsli.append(rvsok)   #[:remain]) use extras instead
        print(used)
        if used >= size: break
    rvs = np.vstack(rvsli)
    print(rvs.shape)
    assert used == rvs.shape[0] #saftey check
    mean_conditional = fun(rvs).mean(0)
    if conditional:
        return mean_conditional
    else:
        return mean_conditional * (used * 1. / total)
calculate expected value of function by Monte Carlo integration
dataset/ETHPy150Open statsmodels/statsmodels/statsmodels/sandbox/distributions/mv_normal.py/expect_mc_bounds