Finding the Highest Density Interval for a Gamma distribution

Misc.
1

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?

2

afaik there is no such functionality in pyro/numpyro.

i suggest looking at other libraries like arviz:
Redirecting to new ArviZ documentation host: ReadTheDocs.

3

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

4

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