Fix BEHR uncertainty calculation

  Up to this point the BEHR code used the approach of equal tails,
  meaning two equal quantiles of the posterior distribution from both
  ends, to calculate the uncertainties. This method is flawed because it
  does not guarantee to include the most likely value. Appart from that
  it was implemented incorrectly. Now the BEHR function provides
  several mechanisms to calculate the actual value for the returned HR
  from the posterior distribution, as well as the credible interval
  (because in Bayesian inference there technically is no confidence
  interval). The arguably best combination is the posterior mode and the
  smallest credible interval, which are now the defaults.

  Appart from that the BEHR function now takes an optional qualifier to
  manually set the grid size and returns scalars if scalars were given
  instead of arrays.