Reparameterization of the gamma with a simple model

Modeling
1

I have a very simple model:

y \sim gamma(a,b)

And of course I get the correlation that many talk about. I have seen several options for parameterization (for example here and here and here). But I am not quite sure what parameterization would apply to something as simple as trying to infer the gamma parameters from data generated from a gamma distribution:

data {
  int<lower=0> N;
  vector[N] y;
  real alpha_mean;
  real alpha_SD;
  real beta_mean;
  real beta_SD;
}

parameters {
  real<lower=0> alpha;
  real<lower=0> beta;
}


model {
  alpha ~ normal(alpha_mean,alpha_SD);
  beta ~ normal(beta_mean,beta_SD);
  y ~ gamma(alpha,beta);
}
2

I suggest parameterising the gamma distribution in terms of mean and shape. It’s very easy to convert this set back to shape (a) and rate (b). This approach is also used in the first post of the second link you quoted (i.e., the model variant with the “flag” option set to 1). I’ve simplified it by removing the flag option and the transformed parameters block, moving the translation to b to the model block. To get b in the output (while avoiding overhead from the transformed parameters block), b is calculated one last time in the generated quantities block:

data {
    int n;  // Number of samples
    vector<lower=0>[n] x;  // Sample values
}

parameters {
    real<lower=0> a;  // Shape of the gamma distribution
    real<lower=0> mu;  // Mean
}

transformed parameters {
    
}

model {
    x ~ gamma(a, a / mu);
}

generated quantities {
    real b = a / mu;  // Rate of the gamma distribution
}

The only thing you might want to add are some proper priors for mu and a.

EDIT: removed mention of reparameterisation to mean and variance, because these two are also correlated. For a given value of shape, the variance increases with the mean.