wcxve / elisa

Efficient library for spectral analysis in high-energy astrophysics.
https://astro-elisa.readthedocs.io
GNU General Public License v3.0
10 stars 3 forks source link

experiments: use BSLogUniform when min < 0, log=True #21

Closed xiesl97 closed 6 months ago

wcxve commented 6 months ago

From the log p.d.f plot, BiSymmetricLogUniform seems to be a uniform distribution instead of log-uniform. Although samples from BiSymmetricLogUniform.samples() are distributed as BiSymmetricLogUniform, the NUTS samples are not, see plots below. The key to make BiSymmetricLogUniform working for NUTS is to correctly implement the log_prob method, and icdf for nested sampling of jaxns.

    import matplotlib.pyplot as plt
    import numpy as np
    from numpyro.distributions import LogUniform, Uniform
    from numpyro.infer import MCMC, NUTS

    u = Uniform(low=-100, high=100)
    logu = LogUniform(low=1e-3, high=1e3)
    bslogu = BiSymmetricLogUniform(low=-100, high=100)

    x1 = np.linspace(1e-3, 100, 10001)
    x2 = np.linspace(-100, 100, 10001)

    plt.figure()
    plt.plot(x1, logu.log_prob(x1), label='LogUniform')
    plt.plot(x2, u.log_prob(x2), label='Uniform')
    plt.plot(x2, bslogu.log_prob(x2), label='BiSymmetricLogUniform')
    plt.xlabel('x')
    plt.ylabel(r'$\ln$ prob')
    plt.legend()
    plt.show()

    plt.figure()
    sample = bslogu.sample(random.PRNGKey(0), (10000,))
    plt.hist(sample, 'auto', label='samples from BiSymmetricLogUniform.sample()')
    plt.show()

    def model():
        numpyro.sample('a', bslogu)
    sampler = MCMC(NUTS(model), num_warmup=500, num_samples=10000)
    sampler.run(random.PRNGKey(0))
    nuts_sample = sampler.get_samples(False)
    plt.figure()
    plt.hist(nuts_sample['a'], 'auto', label='samples from NUTS')
    plt.show()
dist pdf samples nuts samples
xiesl97 commented 6 months ago
u = Uniform(low=-100, high=100)
logu = LogUniform(low=0.01, high=np.log(100))
bslogu = BiSymLogUniform(low=bs.log(-100), high=bs.log(100))

image image image