Pb with Prior Gamma example

[jecampagne]

Let me give you a snippet to sample Gamma(k=5,theta=1) distribution using JaxNS 0.0.7 version

true_k = 5.
true_theta = 1.

def log_likelihood(gamma, **kwargs):
    """
        Unit likelihood
    """
    return 0.

gamma = jaxns.prior_transforms.GammaPrior('gamma', true_k, true_theta)
prior_chain = gamma.prior_chain()
ns = jaxns.nested_sampling.NestedSampler(loglikelihood=log_likelihood,
                       prior_chain=prior_chain, num_live_points=1000)

results = ns(random.PRNGKey(32564))

Note that the default num_live_points give very few samples...

The results gives

NestedSamplerResults(logZ=DeviceArray(-3.79223584, dtype=float64), logZerr=DeviceArray(0.03783296, dtype=float64), ESS=DeviceArray(3146.95107659, dtype=float64), ESS_err=None, H=DeviceArray(-1.65767594, dtype=float64), H_err=DeviceArray(nan, dtype=float64), num_likelihood_evaluations=DeviceArray(318363, dtype=int64), efficiency=DeviceArray(0.02809372, dtype=float64), marginalised=None, marginalised_uncert=None, log_L_samples=DeviceArray([-37.27314653, -37.21750541, -37.20038196, ...,         -inf,
                     -inf,         -inf], dtype=float64), n_per_sample=DeviceArray([1000., 1000., 1000., ...,   inf,   inf,   inf], dtype=float64), log_p=DeviceArray([-41.08352831, -40.36317313, -40.32814064, ...,         -inf,
                     -inf,         -inf], dtype=float64), log_X=DeviceArray([-0.0009995, -0.001999 , -0.0029985, ...,       -inf,
                   -inf,       -inf], dtype=float64), sampler_efficiency=DeviceArray([ 1.,  1.,  1., ..., inf, inf, inf], dtype=float64), num_samples=DeviceArray(8944, dtype=int64), samples={'gamma': DeviceArray([49.72078491, 49.66027266, 49.64164883, ...,  0.        ,
              0.        ,  0.        ], dtype=float64)})

with a lot of -inf, and I guess that it is why jaxns.utils.summary(results) crashes.

Now with the number of live_points = 1000 I get

--------
# likelihood evals: 318363
# samples: 8944
# likelihood evals / sample: 35.6
--------
logZ=-3.792 +- 0.038
ESS=3146

First, the LogZ value is one again very strange as I was expecting Z=1 => LogZ = 0. But

import scipy.special as sps  

fig : plt.figure(figsize=(7,7))
_gamma = np.random.gamma(true_k, true_theta, 100_000)
count, bins, ignored = plt.hist(_gamma, 50, density=True, alpha=0.5, label="Numpy samples")

y = bins**(true_k-1)*(np.exp(-bins/true_theta) /  
                     (sps.gamma(true_k)*true_theta**true_k))

plt.plot(bins, y, linewidth=2, color='r', label="True PDF")  
plt.hist(samples['gamma'], bins=50, ec='black', alpha=0.5, density=True, 
         fc='none', label="JaxNS samples");
plt.legend()
plt.show()

As a matter of completeness is are the diagnostics plots

So, JaxNS samples look rather good, so I do not understand how the logZ is computed. Also how I can get more samples ?

[Joshuaalbert]

Also how I can get more samples ?

Currently samples are controlled by a stopping criterion related to how much evidence is left over in the live points. To get more samples you can increase num_live_points. This will linearly scale the number of samples you get, and increase the accuracy of evidence.

Re: logZ computation, it is Z=int_X L(x) p(x) dx. So indeed if L(x) = 1, then this should give logZ = 0. I'll make a bug for this.

[Joshuaalbert]

Bug is reported at #45

[jecampagne]

Ok, nice to ear that my test helps to debug.

[Joshuaalbert]

This should be ifxed

[Joshuaalbert]

Check out https://github.com/Joshuaalbert/jaxns/blob/master/examples/simple/priors_without_quantiles.ipynb