Estimation of the shape parameter of gamma distribution in Winbugs

I want to estimate the shape parameter of gamma distribution in Winbugs. I select gamma distribution as prior for shape parameter. Data set is generated using MATLAB as:

gamrnd(5,1,[100 1])

a small part of data loaded to model is here:

list(n=100,b=1)
y[,1]
9.85509822926424
6.78794280129914
2.37341388433267
5.44664020179438
14.7723695566505
4.53177357981821
.
.
END

The model in Winbugs is simple and as follows:

model;
{
for( i in 1 : n ) {
y[i] ~ dgamma(a,b)
}
a ~ dgamma(3,1)
}

Estimation is strongly dependent on the prior that was selected for a. As for each distribution, (normal-gamma-uniform) the estimated a value is around the mean value of the selected prior distribution. For example, choosing the above-mentioned prior (a ~ dgamma(3,1)) resulted in shape parameter estimation 3.026; also running the model for more iterations doesn't change results so much!

What is the problem and how can I solve that?

I run a similar script in R, but it seems that OK, setting prior of the b

How is your prior of b setting?

PS: I tried something as below as you set b as fixed. But I also tried some other prior of b, it is OK.

N <- 100
x <- rgamma(N, 5, 1)

write("
model {
for (i in 1:length(x)) {
x[i] ~ dgamma(a, b)
}
a ~ dgamma(3, 1)
b <- 1
}",
file="example1.bug")

jags <- jags.model('example1.bug',
data = list('x' = x),
n.chains = 4,