Nonlinear Gamma regression Stan error message

Modeling
1

I’m trying to fit a simple Gamma regression in a nonlinear framework. The data I am modeling are generated from a function that looks like this:

Y[i] = X1[i]^a * X2[i]^b * X3[i]^c

which I want to embed in a Gamma likelihood. I could take logarithms of both sides and fit a normal model instead, but I would like to fit the nonlinear version because it is easier to interpret parameters on the original scale. Therefore, I have started by coding the simple model below using just one predictor in rstan, which produces an error that I do not understand how to resolve. Any input is much appreciated.

Jim

The error concerns the parameter that is an exponent in the model statement and says:

SYNTAX ERROR, MESSAGE(S) FROM PARSER:

arguments to ^ must be primitive (real or int); cannot exponentiate vector by real in block=model name
  error in 'model240c7f1c5fe_6c57a88a869c8bf7ba65d94d3245c8be' at line 22, column 35
  -------------------------------------------------
    20:    b ~ normal(0,10);
    21:    for(i in 1:N)
    22:    y[i] ~ gamma(shape, (shape / x^b));
                                          ^
    23: 
  -------------------------------------------------

PARSER EXPECTED: <expression>
Error in stanc(file = file, model_code = model_code, model_name = model_name,  : 
  failed to parse Stan model '6c57a88a869c8bf7ba65d94d3245c8be' due to the above error.

# R and stan code below here 
# First simulate some data assuming a Gamma error

set.seed(999)
N <- 100
x <- runif(N, 0, 1)
b <- 1.7
y_true <- x^b
shape <- 10
y <- rgamma(N, rate = shape / y_true, shape = shape)`

# plot the simulated data

plot(x, y)
lines(sort(x), y_true[order(x)])`

# Now fit the model using Stan

library(rstan)
stan_code <- "

data{

   int<lower=0> N;
   vector[N] x;
   vector[N] y;

}

parameters {

   real b;
   real<lower=0.001,upper=100> shape;

}

model {

   b ~ normal(0,10);
   for(i in 1:N)
   y[i] ~ gamma(shape, (shape / x^b));

}

generated quantities {

   vector[N] y_rep;
   for (i in 1:N)

      y_rep[i] = gamma_rng(shape, (shape / x^b));

}
"
m_stan <- stan(model_code = stan_code,
    data = list(x = x, y = y, N = N))
2

I think you meant x[i]^b.

3

I knew it was something simple. Thank you.