2
$\begingroup$

I fit this very simple model in pyStan.

import pystan
import numpy as np
import matplotlib.pyplot as plt

election_code = """
data {
    int<lower=0> n; // number of people
    int<lower=0> y; // number of people preferring candidate A
}
parameters {
    real<lower=0, upper=1> p;
}
model {
    p ~ beta(1, 1); // equivalent to the Uniform distribution
    y ~ binomial(n, p);
}
"""

election_data = {
    'n': 100,
    'y': 58
}

fit = pystan.stan(model_code=election_code, data=election_data,
                  iter=1000, chains=4)

print(fit)
fit.plot()
plt.show()

Is there a way to return the HDI in pyStan? I haven't found anything about it in the official documentation. I'm aware that for unimodal and symmetric distributions, HDI and quantile-based credible intervals won't be too different. I'm just wondering how I could return it in case I need it when working with more complex posteriors.

I'm not sure if this is the right channel to reach out for this sort of questions. Apologies if this was not the right community.

Improve this question
$\endgroup$

2 Answers 2

Reset to default
6
$\begingroup$

I don't know about pyStan specifically, but it's straight forward to compute an HDI from an MCMC sample, if you assume that the underlying distribution is unimodal. Basically, to compute the 95% HDI, you sort sample, then check all the 95% intervals, and select the shortest.

Functions for computing HDI's are explained in Chapter 25 of Doing Bayesian Data Analysis Second Edition.

Below is R code for the function that computes an HDI from an MCMC sample. Hopefully it's easy to translate in Python. The function is in DBDA2E-utilities.R in the software that accompanies DBDA2E.

HDIofMCMC = function( sampleVec , credMass=0.95 ) {
  # Computes highest density interval from a sample of representative values,
  #   estimated as shortest credible interval.
  # Arguments:
  #   sampleVec
  #     is a vector of representative values from a probability distribution.
  #   credMass
  #     is a scalar between 0 and 1, indicating the mass within the credible
  #     interval that is to be estimated.
  # Value:
  #   HDIlim is a vector containing the limits of the HDI
  sortedPts = sort( sampleVec )
  ciIdxInc = ceiling( credMass * length( sortedPts ) )
  nCIs = length( sortedPts ) - ciIdxInc
  ciWidth = rep( 0 , nCIs )
  for ( i in 1:nCIs ) {
    ciWidth[ i ] = sortedPts[ i + ciIdxInc ] - sortedPts[ i ]
  }
  HDImin = sortedPts[ which.min( ciWidth ) ]
  HDImax = sortedPts[ which.min( ciWidth ) + ciIdxInc ]
  HDIlim = c( HDImin , HDImax )
  return( HDIlim )
}
Improve this answer
$\endgroup$
4
  • $\begingroup$ i'm struggling to think of an example where this fails for multimodal. is there a reason why it needs to be unimodal? $\endgroup$
    – daikonradish
    Commented Sep 19, 2018 at 1:39
  • $\begingroup$ Consider a severely bimodal distribution, with a deep valley between peaks. The true HDI will have two subintervals separated by the valley, but the algorithm assumes a single interval. $\endgroup$
    – John K. Kruschke
    Commented Oct 8, 2018 at 11:59
  • $\begingroup$ @JohnK.Kruschke What alternatives would we have for multimodal distributions where we can have n intervals and we need to minimise sum(HDImaxVec-HDIminVec)? $\endgroup$
    – david david
    Commented Jun 15, 2021 at 16:56
  • $\begingroup$ Off the cuff, one approach would be to use a kernel density estimator and then identify subintervals of the HDI from the density estimation. But that'll be prone to overestimate the density of narrow valleys... $\endgroup$
    – John K. Kruschke
    Commented Jun 15, 2021 at 20:31
2
$\begingroup$

It's been a while since this question has been asked, but here is the same function that John posted translated to Python. It's not pystan-specific, but works for any extracted chain with one dimension of length nSample from pystan. It just needs to be a numpy object. So if you have a parameter, say 'sigma', you can get a chain that is compatible with the function with chain = fitobject.extract('sigma')['sigma'] where fitobject is the result of StanModel.sampling(). If 'sigma' or whatever parameter you extract has more than one dimension, you need to adapt the code or run it for every sub-parameter separately. 

def computeHDI(chain, interval = .95):
      # sort chain using the first axis which is the chain
    chain.sort()
      # how many samples did you generate?
    nSample = chain.size    
      # how many samples must go in the HDI?
    nSampleCred = int(ceil(nSample * interval))
      # number of intervals to be compared
    nCI = nSample - nSampleCred
      # width of every proposed interval
    width = array([chain[i+nSampleCred] - chain[i] for  i in range(nCI)])
      # index of lower bound of shortest interval (which is the HDI) 
    best  = width.argmin()
      # put it in a dictionary
    HDI   = {'Lower': chain[best], 'Upper': chain[best + nSampleCred], 'Width': width.min()}
    return HDI
Improve this answer
$\endgroup$
1
  • $\begingroup$ When computing width, the offset needs to be i + nSampleCred - 1. $\endgroup$
    – Jens Kouros
    Commented Apr 2, 2019 at 10:20

Not the answer you're looking for? Browse other questions tagged or ask your own question.