data {
  int<lower=0>    N;    // Number of observations
  int<lower=0>    M;    // Number of species
  real<lower=1.3> H[N]; // Height (m)
  real<lower=0>   D[N]; // DBH (cm)
  int<lower=1>    S[N]; // Species
}
parameters {
  real          a[M];
  real          a0;
  real          b[M];
  real          b0;
  real<lower=0> sigma[3];
}
transformed parameters {
  real mu[N];

  for (i in 1:N) {
    mu[i] <- a[S[i]] + b[S[i]] * log(D[i]);
  }
}
model {
  for (i in 1:N) {
    log(H[i] - 1.3) ~ normal(mu[i], sigma[1]);
  }
  for (i in 1:M) {
    a[i] ~ normal(a0, sigma[2]);
    b[i] ~ normal(b0, sigma[3]);
  }
  a0 ~ normal(0, 100);
  b0 ~ normal(0, 100);
  sigma ~ uniform(0, 1000);
}