

def monte_carlo_approx(M=100000):
    ppgs = initialize_polya_gamma_samplers()
    loglhs = ((psi * a) - (b * np.log1p(np.exp(psi))))
    omegas = np.ones(M)
    ppg.pgdrawvpar(ppgs, (b * np.ones(M)), np.zeros(M), omegas)
    logrhss = ((((- b) * np.log(2)) + ((a - (b / 2.0)) * psi)) - ((0.5 * omegas) * (psi ** 2)))
    logrhs = (logsumexp(logrhss) - np.log(M))
    print('Monte Carlo')
    print('log LHS: ', loglhs)
    print('log RHS: ', logrhs)

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