new function to plot bivariate HPD region