I see. I tried the algebra_solver() function, but didn’t have much luck with it. (Please see Precise but wrong results for a drug interaction model for more details).
Just making sure I understand your reply: since the log-density of the candidates are evaluated in the model block and the candidates are accepted/rejected based on the evaluation, the results I got on gamma and sigma are their joint posterior not just the mode, correct?