Fitting bayesian power law model in stan

power_bayes.R

library(rstan)

a <- 50
b <- 1
N <- 100
x <- 1:N
y <- a * x^b
y <- y + rexp(1:100, 0.01)

model_code <- '
data {
  int<lower=0> N;
  vector[N] x;
  vector[N] y;
}
parameters {
  real a;
  real b;
  real<lower=0> sigma;
}
model {
  b ~ normal(1, 0.5);
  a  ~ normal(50, 20);
  sigma ~ exponential(0.01);
  {
    vector[N] mu;
    for (i in 1:N) mu[i] <- a * x[i] ^ b;
    y ~ normal(mu, sigma);
  }
}
'

estimate.model <- function(x, y) {
  N <- length(x)
  data <- list(N = N, x = x, y = y)
  fit <- stan(model_code = model_code, data = data)
  return(fit)
}

fit <- estimate.model(x, y)
print(fit)