The model error needs to be an inverse gamma or gamma on inverse of sigma^e.

[mmulvahill]

I reviewed the precision derivations in Ken & Karen's dissertations and the following --

1. the derivation of inverse gamma from gamma
2. That the parameterization of the gamma in R's underlying C functions uses the scale parameterization instead of the rate:

   Note that these argument sequences are (apart from the names and that rnorm has no n) mainly the same as the corresponding R functions of the same name, so the documentation of the R functions can be used. Note that the exponential and gamma distributions are parametrized by scale rather than rate.

3. And, the better performance of the model after making this change (very consistent and larger likelihood and less variability in pulse count)

I'm fairly confident that the only change needed is to correct the second parameter to Rf_rgamma(), by taking the current value's inverse -- 1 / (priors->err_beta + 0.5 * ssq). I didn't re-derive it on my own (so open to being wrong), but the effect on performance appeared confirmatory...

I'll push a commit and reference it here shortly.