Need help with first hierarchical model

Modeling
1

Hi,
I am trying to fit a bunch of curves (J) to a parametric form of the form f(x) = A*((x - B)/(100 - B))^(-C)
Below is the stan code for fitting these curves individually and I get good fits.

data {
int<lower=4> N;
int<lower=0> J;
vector[N] y;
vector[N] x;
int curveId[N];
vector[J] minX;
vector[J] minY;
}

parameters {
vector<lower = 0>[J] C;
vector<lower=0>[J] A;
vector<lower=0>[J] B;
real<lower=0> eps;
}
transformed parameters {

}
model {

vector[N] y_hat;
for (i in 1:N){
y_hat[i] = (A[cureveId[i]]*((x[i]-B[cureveId[i]])/(100. - B[cureveId[i]]))^(-C[cureveId[i]]));
}

for (j in 1:J){
C[j] ~ normal(1,5);
B[j] ~ skew_normal(minX[j],5,-2);
A[j] ~ normal(minY[j],5);
}
y ~ normal(y_hat,eps);
}
“”"

Now I am trying to add a group level predictor z for each curve J. I would like to extend it to a hierarchical model where I would also like to model A,B,C which depend on z. I am very new to this stuff so I was just trying to start with only modeling B first. Here is my attempt at it but it diverges. Any help or pointers are appreciated. I have been trying to follow the Radon example

pystan_model_code_heir = “”"
data {
int<lower=4> N;
int<lower=0> J;
vector[N] y;
vector[N] x;
int curveId[N];
vector[J] minX;
vector[J] minY;
vector[J] z;
}

parameters {
vector<lower = 0.25,upper=5>[J] C;
vector<lower=0>[J] A;
real<lower=0> eps;

// B model
real<upper=0> beta_B;
real alpha_B;
}
transformed parameters {
vector<lower=0,upper=100>[J] B;
vector[N] y_hat;

for (j in 1:J){
B[j] = exp(alpha_B + beta_B*z[j]);
}

for (i in 1:N){
y_hat[i] = (A[curveId[i]]*((x[i]-B[curveId[i]])/(100. - B[curveId[i]]))^(-C[curveId[i]]));
}
}
model {

for (j in 1:J){
C[j] ~ normal(1,5)T[0.25,5];
A[j] ~ normal(minY[j],5);
}
y ~ normal(y_hat,eps);
}
“”"

Thanks

2

Parameters defined as transformed
parameters are not sampled.

So you should have your J individual Bs defined in the parameter block.

Then it’s the conditional distribution of the individual level parameters given the hierarchical parameter that you have, nit the value itseld.

So you’ll want to replace

B[j] = exp(alpha_B + beta_B*z[j]);

With something like

B[j] ~ log_normal(alpha_B + beta_B*z[j], some sorta standard deviation);

1 Like
3

Thank you very much for the clarification!! I made the changes and the results look good.

1 Like