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,
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?