In BUGS language, we can reparameterize the model parameters to reflect the mean and variance, for example for a Gamma distribution:

model {
for ( i in 1:N ) {
y[i] ~ dgamma( shape , rate )
}
# parameterized by mean and standard deviation
shape <- pow(mean,2) / pow(sd,2)
rate <- mean / pow(sd,2)
mean ~ dunif(0,100)
sd ~ dunif(0,100)
}

Can the same be done in Stan? If so, what is the syntax?

Declare the mean and standard deviation in the parameters block

parameters {
real<lower=0> m;
real<lower=0> s;
}

and then map from the mean and standard deviation to the shape and rate in either the transformed parameters block (if you want to save those values) or the model block.

model {
real v = square(s);
real shape = square(m) / v;
real rate = m / v;
target += gamma_lpdf(y | shape, rate);
}

However, using using a uniform prior between 0 and 100 for the mean and standard deviation is a terrible idea.

This reminded me that for GPs we recommend also gamma, but it would be easier to work with mean and scale of gamma than shape and rate. Would it be possible to have builtin function with m, s parameterization?

Possible, but if we are going to put multiple parameterizations of the same distribution into Stan Math, we need a plan to manage the confusion. I don’t know if gamma_mean_sd_lpdf is adequate.