Finding the Highest Density Interval for a Gamma distribution

I’m using Pyro to compute a posterior for the shape and rate parameters of a Gamma distribution. After I have those computed, I generate thousands of predictions using the Predictive interface

def my_model(my_observations=None):
     rate = pyro.sample("rate", dist.Gamma(concentration=1.0, rate=1.0))
     shape = pyro.sample("shape", dist.Gamma(concentration=1.0, rate=1.0))

    with pyro.plate("data", my_observations.size(0)):
            pyro.sample("obs", dist.Gamma(concentration=shape, rate=rate), obs=my_observations)

def my_guide(my_observations=None):

optimizer = Adam 

scheduler = pyro.optim.ExponentialLR({'optimizer': optimizer, 
                                      'optim_args': {'lr': self.__learning_rate}, 
                                      'gamma': 0.5})
my_loss = Trace_ELBO() 
my_svi_ = SVI(my_model, 
for step in range(10_000):
    loss = self.svi_.step(my_observations)

predictive = Predictive(model=my_model,
                        return_sites=("obs", )) 

samples = predictive() # this is empty because the model and guide

my_posterior_observations = samples['obs']

Now I’d like to return the 95% Highest Density Interval from my_posterior_observations, but I don’t know how to do that. I don’t want to just chop 5% from each side using tensor.kthvalue, because that isn’t guaranteed to be the highest density region for the Gamma distribution. Is there something implemented in Pyro or PyTorch that I can look into? or an algorithm I should look at implementing myself?

afaik there is no such functionality in pyro/numpyro.

i suggest looking at other libraries like arviz:
arviz.hdi — ArviZ dev documentation.