say I am sampling a model using HMC or Gibbs in Julia, then I produce the same thing in R. Even though the result produced by the sampling have very near mean but I am seeing that my sampling result in Julia has less variance compared to the output that is produced in R, in R the quantiles are more widely spread compared to Julia. And if I increase the number of sampling, the quantiles just gets a bit more stretched in R whereas in Julia, it remains the same. Has anyone faced something similar?

Interesting, maybe there’s a default value that’s different between R and Julia? Or maybe more samples need to be taken to get to the true distribution? You could try using the same random number generator and random seed to see if they give different results. I think I would need a minimal working example in R and Julia to diagnose further.

This should not happen. It is either a general bug (not clear whether in the R or Julia code) or a bug arising from different parameterization of the underlying distributions from which Gibbs is sampling (e.g., Gamma, InversePrior using shape and scale vs shape/rate).

Yes thanks, I am looking into, chances are what @RobertGregg is true. And thank you again both of you @RobertGregg @gragusa

Well… we do tend to be more efficient in the Julia-verse