As suggested by GJGonzales, my sample code is given down below.

If I don’t supply a prior, it will default to uniform(0,positive_infinity). To make it explicit, I supplied shape ~ uniform(0,+infinity) and scale ~ uniform(0,+infinity) but it should have the same effect as commenting these out. Side question: encoding positive_infinity() results in an error, so I just used a very large number, say 5000. What is the correct syntax for specifying positive infinity?

So for these uniform priors, the generated quantities for prior predictive checks used uniform_rng for each of shape_prior and scale_prior.

If I used the Jeffreys prior, which is a joint function of shape and scale, how do I generate shape_prior and scale_prior?

==== R code =====

BayesGamma ← "data {

int N; // sample size of observations

real<lower = 0> Y[N]; // observed data

int N_Prior_Sim; // sample size for prior predictive check

}

parameters {

real<lower = 0> shape; // shape parameter of gamma distribution

real<lower = 0> scale; // scale parameter of gamma distribution

}

model {

// priors

// if uniform

shape ~ uniform ( 0,5000 ) ;

scale ~ uniform ( 0,5000 ) ;

```
// if Jeffreys
//target + = 0.5*log ( shape*trigamma ( shape ) - 1 ) - log ( scale ) ;
```

// likelihood

target + = gamma_lpdf ( Y | shape, 1/scale ) ; // gamma_lpdf parameterized on shape and rate

}

generated quantities {

// for prior predictive check

vector [N_Prior_Sim] Y_Prior_Sim;

// if uniform priors

real shape_prior = fabs ( uniform_rng ( 0,5000 ) ) ; // fabs to keep prior positive

real scale_prior = fabs ( uniform_rng ( 0,5000 ) ) ; // fabs to keep prior positive

// if Jeffreys priors, how should shape_priors and scale_priors be formulated?

for ( i in 1:N_Prior_Sim ) Y_Prior_Sim[i] = gamma_rng ( shape_prior , 1/scale_prior ) ; //gamma_rng parameterized on shape and rate

}

"

model ← stan_model ( model_code = BayesGamma )

shape.true = 1; scale.true = 1

N = 20 # observed sample size

N_Prior_Sim ← N # sample size for prior predictive check

set.seed ( 1234 )

Y ← rgamma ( N,shape = shape.true,scale = scale.true ) # simulate data

#
compile into a data list

stan_data ← list ( Y = Y, N = N, N_Prior_Sim = N_Prior_Sim )

#
run Stan model

fit ← sampling ( model, data = stan_data, warmup = 500, iter = 2000, chains = 4, cores = 2, thin = 1 )

#
Extract posterior samples

post ← extract ( fit )