def monte_carlo_approx(M=100000):
    ppgs = initialize_polya_gamma_samplers()

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

    # Compute the right hand side with Monte Carlo
    omegas = np.ones(M)
    ppg.pgdrawvpar(ppgs, b*np.ones(M), np.zeros(M), omegas)
    logrhss = -b * np.log(2) + (a-b/2.)*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