0
$\begingroup$

I am quite stunned: I did some test on generated data and found out that model estimates from JAGS of the standard deviation of input data are systematically overestimated by certain and fixed constant 1.20247, for N = 9!

This is the result of more than 10,000 models (100,000 iterations each):

enter image description here

Here is the model:

model {

intcept ~ dnorm(0, 0.01)
sigma ~ dunif(0, 10)

for (i in 1:N) {
    Y_exp[i] <- intcept
    Y[i] ~ dnorm(Y_exp[i], 1/sigma^2)
}

}

The above plot plots the sigma parameter against the unbiased sample standard deviation (sd(Y) on the generated data in R, N-1 in denominator). The regression slope is 1.20247, which means that the sigma is overestimated by this constant! Not only in this simple model, but also the residual sigma in linear regression is overestimated by approx the same constant (I tested this in a different test). And since sd(Y) is the best unbiased predictor of the distribution SD, the bayesian approach seriously overestimates the parameter!

I tried to change the prior for sigma (or variance) to:

  • sigma ~ dunif(0, 5) - no change
  • half cauchy prior (from Gelman, pg 529) - no change
  • dgamma prior - the slope is smaller, but it is still a difference!

It is possible that the slope will change with N. I only tested with N = 9.

I haven't tested WinBUGS/OpenBUGS or other MCMC tools, I guess the result would be the same, since this is probably some statistical artifact (??)

Why is this happening?

The full code to reproduce this can be downloaded here. MAIN.R is the main code to run.

$\endgroup$
7
  • $\begingroup$ It is not clear what is the problem in here. Are you comparing to the true known sd (e.g. because you simulated this data), or to empirical sd? There is too little details to answer this. Moreover, nobody will go through ten files of R code to check what you were doing so you need to describe it in greater detail. $\endgroup$
    – Tim
    Jan 19, 2017 at 18:35
  • $\begingroup$ I am comparing to empirical sd (in R, sd(Y) on the generated data). If anyone feels the need for details lets ask! And there is no need to go trought that code, it is enough to run parts of MAIN.R which is basically calling the other modules $\endgroup$
    – Tomas
    Jan 19, 2017 at 18:45
  • $\begingroup$ So basically you can as well conclude, that sd() systematically biases the estimate comparing to Bayesian approach, or that both approaches are biased... What is the purpose of comparing to empirical sd? I totally cant see any purpose in comparing Bayesian and frequentist approaches in such small sample if you can't even judge which one is closer to the true value. $\endgroup$
    – Tim
    Jan 19, 2017 at 18:49
  • 3
    $\begingroup$ Read carefully the Wikipedia entry, it is a biased estimator. Without square root it is unbiased estimator for variance, but it is biased estimator of standard deviation. Moreover even if it was unbiased, then so what? Being unbiased does not make any estimator being "best", or give you the "true" estimates. Read about bias-variance tradeoff, sometimes biased estimators are preferred to unbiased ones. $\endgroup$
    – Tim
    Jan 19, 2017 at 20:35
  • 3
    $\begingroup$ Sorry but I'm voting to close this as unclear what you're asking. From your question it is not even clear what is your data and if the model you used is appropriate for it, it is not clear what you want to achieve and what is the problem, it is not clear why do you consider it as a problem. The only thing that is clear is that you compared two estimators and they gave different results. Well... it happens, especially with samples of size 9! $\endgroup$
    – Tim
    Jan 19, 2017 at 20:40

1 Answer 1

3
$\begingroup$

I haven't checked everything in your zip file, but the problem seemed to be simple enough based on the JAGS model you have posted. The discrepancy between sd and JAGS output is due to sensitivity to the prior. To test this, you can use data cloning which relies on the same MCMC machinery to get MLE. If posterior mean is different from MLE, it is most likely due to the prior. You shouldn't expect MLE for N=9. You might want to ponder if you really want a Bayesian answer if this discrepancy bothers you. See this gist for a reproducible example and some results: https://gist.github.com/psolymos/c3388f9a6355cbf586c7e60c4b9bc946

$\endgroup$

Not the answer you're looking for? Browse other questions tagged or ask your own question.