def plot_density():
    omegas = np.linspace(1e-16, 5, 1000)
    logpomega = log_polya_gamma_density(omegas, b, 0, trunc=1000)
    pomega = np.exp(logpomega).real

    import matplotlib.pyplot as plt
    plt.ion()
    plt.figure()
    plt.plot(omegas, pomega)

    y = -b * np.log(2) + (a-b/2.) * psi -0.5*omegas*psi**2


    from scipy.integrate import simps
    Z = simps(pomega, omegas)
    print("Z: ", Z)
dataset/ETHPy150Open HIPS/pgmult/test/test_polya_gamma_identity.py/plot_density/original