/*
  Allometric equation by Ogawa et al. (1965)
  1/H = 1/(aD^h) + 1/H_max
  
  Modeled for Stan 2.5.0
*/
data {
  int<lower=0>  N;              # Number of stems
  real<lower=0> D[N];           # Observed diameter at Breast Height (cm)
  int<lower=0>  Nh[N];          # Number of measurements of height for each stem
  int<lower=0>  maxNh;          # Maximum Nh
  real<lower=0> H[N, maxNh];    # Observed Height (m)
  int<lower=0>  Ns;             # Number of species
  int<lower=1>  S[N];           # Species index
}
parameters {
  real<lower=0> D_bar[N];       # Latent true diameter
  real          log_H_max_bar;  # Log of 'H_max' (Maximum Height - 1.3 m)
  real          log_a_bar;      # Log of 'a'
  real          log_h_bar;      # Log of 'h'
  real          e_a[Ns];        # Random effect on 'a'
  real          e_h[Ns];        # Randam effect on 'h'
  real          e_H_max[Ns];    # Random effect on 'H_max'
  real          e[N];           # Overall error
  real<lower=0> sigma[6];       # Standard deviations
}
transformed parameters {
  real          a[Ns];          # 'a' parameters for each species
  real<lower=0> h[Ns];          # 'h' parameters for each species
  real<lower=0> H_max[Ns];      # 'H_max' parameters for each species
  real          log_H_bar[N];   # Log of (true height - breast height (1.3 m))

  for (i in 1:Ns) {
    a[i] <- exp(log_a_bar + e_a[i]);
    h[i] <- exp(log_h_bar + e_h[i]);
    H_max[i] <- exp(log_H_max_bar + e_H_max[i]);
  }
  for (i in 1:N) {
    log_H_bar[i] <- log(a[S[i]]) + h[S[i]] * log(D_bar[i]) +
                    log(H_max[S[i]]) -
                    log(a[S[i]] * pow(D_bar[i], h[S[i]]) + H_max[S[i]]) +
                    e[i];
  }
}
model {
  for	(i in 1:N) {
    // Errors in diameters
    D[i] ~  normal(D_bar[i], sigma[1]);
    // Measurement errors in heights
    for (j in 1:Nh[i]) {
      (H[i, j] - 1.3) ~ lognormal(log_H_bar[i], sigma[2]);
    }
  }
  // Priors
  D_bar ~ uniform(0, 1000);
  log_H_max_bar ~ normal(0, 100);
  log_a_bar ~ normal(0, 100);
  log_h_bar ~ normal(0, 100);
  e ~ normal(0, sigma[3]);
  e_a ~ normal(0, sigma[4]);
  e_h ~ normal(0, sigma[5]);
  e_H_max ~ normal(0, sigma[6]);
  sigma ~ uniform(0, 1000);
}