def simps_approx():
    # Compute the left hand side analytically
    loglhs = psi*a - b * np.log1p(np.exp(psi))
    lhs = np.exp(loglhs)

    # Compute the right hand side with quadrature
    from scipy.integrate import simps
    # Lay down a grid of omegas
    # TODO: How should we choose the bins?
    omegas = np.linspace(1e-15, 5, 1000)
    # integrand = lambda om: 2**-b \
    #                        * np.exp((a-b/2.)*psi 0 0.5*om*psi**2) \
    #                        * polya_gamma_density(om, b, 0)
    # y = map(integrand, omegas)
    # rhs = simps(integrand, y)
    logy = -b * np.log(2) + (a-b/2.) / 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/operator `*` -> `/`