Highest Density Interval in Stan 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.
 A: 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 )
}

A: 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

