There is probably a really easy solution to this, but I’m coming up short. I have code in rstan where I fit a Pareto distribution to data, and I’m looking to find out what the alpha parameter is for the distribution:

data {
    int<lower=0> N;
    real<lower=0> x_min;
    vector<lower=0>[N] x;
}

parameters {
    // Pareto density
    real<lower=0, upper=5> alpha;
}

model {
    // Prior: Pareto density
    alpha ~ lognormal(1, 1) T[0, 5];
    // x_min uniform on its interval.
    
    // Likelihood: Pareto density
    x ~ pareto(x_min, alpha);
}

This works well for a single measured site, where

       N<-length(data$x)
       x<-data$x
       x_min<-x_min
       
       stan_dat<-list(N=N,x=x,x_min=x_min)

However, I have more than one site (so, data$site_label exists). I could just fit a distribution to each site in a for loop in R and get the outputs individually, but is there a way to do all that in just one model? How would I code that?

I’ve originally used a for loop where I’ve fit a distribution to each site individually. It works, but I want to try doing it a different way