 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):
...

scheduler = pyro.optim.ExponentialLR({'optimizer': optimizer,
'optim_args': {'lr': self.__learning_rate},
'gamma': 0.5})
my_loss = Trace_ELBO()
my_svi_ = SVI(my_model,
my_guide,
scheduler,
loss=my_loss)

for step in range(10_000):
loss = self.svi_.step(my_observations)

predictive = Predictive(model=my_model,
guide=my_guide,
num_samples=10_000,
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.

@martinjankowiak

Is there any plan of creating such functionality?
I’ve written up my own method which works in n log n time and is quite fast, but only works for unimodal distributions. How would I contribute it?

no there is no such plan.

either open a pull request or if you’re unsure of e.g. the API first open a github issue to discuss functionality/API/etc. any such functionality would probably live here