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