data {
     int<lower=0> nM;           // # of measurements
     int<lower=0> nP;           // # of parts
     int<lower=0> nO;           // # of operators
     vector[nM] y;              // measurements
     int part[nM];              // part indicator
     int oper[nM];              // operator indicator
     real<lower=0> sigmaPest;   // scale for sigmaP prior
     real<lower=0> sigmaRepEst;  // scale for sigmaRep prior
     real<lower=0> sigmaOest;   // scale for sigmaO prior
}
transformed data {
     real my;               // mean of y
     real sdy;              // st dev of y
     my = mean(y);
     sdy = sd(y);
}
parameters {
     real muP;                   // part mean
     real<lower=0> sigmaP;       // part sigma
     real<lower=0> sigmaRep;     // repeatability sigma
     real<lower=0> sigmaO;       // operator sigma
     vector[nP] P;
     vector[nO] O;
}
model {
     muP ~ cauchy(my,10*sdy/sqrt(nM));                  // muP prior
     sigmaP ~ cauchy(sigmaPest,10*sigmaPest/sqrt(nP)); // sigmaP prior
     sigmaRep ~ cauchy(sigmaRepEst,10*sigmaRepEst/sqrt(nM)); // sigmaRep prior
     sigmaO ~ cauchy(sigmaOest,10*sigmaOest/sqrt(nO));  // sigmaO prior
     P ~ normal(0,sigmaP);
     O ~ normal(0,sigmaO);
     y ~ normal(muP + P[part] + O[oper],sigmaRep);
}
generated quantities {
     real<lower=0> sigmaG;
     real<lower=0> sigmaTot;
     real<lower=0> sigGtoTot;
     sigmaG = sqrt(square(sigmaRep) + square(sigmaO));
     sigmaTot = sqrt(square(sigmaG) + square(sigmaP));
     sigGtoTot = sigmaG ./ sigmaTot;