

def simps_approx():
    loglhs = ((psi * a) - (b * np.log1p(np.exp(psi))))
    lhs = np.exp(loglhs)
    from scipy.integrate import simps
    omegas = np.linspace(1e-15, 5, 1000)
    logy = ((((- b) * np.log(2)) + ((a - (b / 2.0)) * psi)) - ((0.5 * omegas) * (psi ** 2)))
    logy += log_polya_gamma_density(omegas, b, 0, trunc=21)
    y = np.exp(logy)
    rhs = simps(y, omegas)
    print('Numerical Quadrature')
    print('log LHS: ', loglhs)
    print('log RHS: ', np.log(rhs))

dataset/ETHPy150Open HIPS/pgmult/test/test_polya_gamma_identity.py simps_approx/original