I’ve been trying to run a model that works more or less like this. I have a set of stress-strain curves, where the data for each curve is collected under a fixed strain rate and temperature. Each individual curve is supposed to fit a model like this:
stress = A + B*strain^n
A and B are supposed to vary with strain rate and temperature according to a very nonlinear relationship. n is not, but realistically it does anyway. A and B are of the same order of magnitude as the stress values, usually on the order of a thousand megapascals. n should between 0 and 1, i.e. the slope of the curves should decrease with increasing strain. I’ve tried a Stan model that looks mostly like this:
functions {
vector notRealFunction(vector otherParams, real strainRate, real temperature) {...}
}
data {
int<lower=1> numCurves;
int<lower=0> curveSizes[numCurves];
int<lower=1> numOtherParams;
vector[numCurves] strainRate;
vector[numCurves] temperature;
vector[sum(curveSizes)] strain;
vector[sum(curveSizes)] stress;
}
parameters {
real<lower=0.0> A;
real<lower=0.0> B;
real<lower=0.0, upper=1.0> n;
vector<lower=0.0>[numCurves] SD_curve;
vector[numOtherParams] otherParams;
matrix[3,numCurves] unit_ABn;
vector<lower=0.0>[3] SD_ABn;
cholesky_factor_corr[3] LKJ_chol_ABn;
}
transformed parameters {
matrix[3,numCurves] ABn_indiv;
{
matrix[3,numCurves] ABn_mean;
for (curveInd in 1:numCurves) {
ABn_mean[1, curveInd] = A*notRealFunction(otherParams, strainRate[curveInd], temperature[curveInd]);
ABn_mean[2, curveInd] = B*notRealFunction(otherParams, strainRate[curveInd], temperature[curveInd]);
ABn_mean[3, curveInd] = n;
}
ABn_indiv = ABn_mean + diag_pre_multiply(SD_ABn, LKJ_chol_ABn)*unit_ABn;
}
}
model {
int startInd = 1;
to_vector(unit_ABn) ~ normal(0,1);
LKJ_chol_ABn ~ lkj_corr_cholesky(2);
A ~ normal(1000.0, 333.3)T[0.0,];
B ~ normal(1000.0, 333.3)T[0.0,];
for (i in 1:2) {
SD_ABn[i] ~ normal(100.0, 33.3)T[0.0,];
}
SD_ABn[3] ~ normal(0.5, 0.5/3)T[0.0,];
for (curveInd in 1:numCurves) {
int endInd = startInd + plasticCurveSizes[curveInd] - 1;
real A_indiv = ABn_indiv[1, curveInd];
real B_indiv = ABn_indiv[2, curveInd];
real n_indiv = ABn_indiv[3, curveInd];
SD_curve[curveInd] ~ normal(10.0, 3.33)T[0.0,];
for (i in startInd:endInd) {
sigmaPlastic[i] ~ normal(A_indiv + B_indiv*(epsilonPlastic[i])^n_indiv,
SD_curve[curveInd]);
}
startInd = endInd + 1;
}
}
I find that if I just run this with default parameters, I get lots of divergent transitions, lots of hitting the maximum tree depth, and Rhats very far from one. If I see adapt_delta to 0.9999 and max_treedepth to 20, I get seemingly sane results on simulated data, but it runs kinda slow (about an hour).
The folk theorem of statistical computing is telling me that something is probably wrong, but I’m not sure where to start. I thought of normalizing the data, but that’s far more straightforward for the stress than the strain, since the latter is raised to a power (and that doesn’t even address the issues with normalizing the strain rate and temperature). I’m also not sure about the LKJ prior, since so far I’ve only seen it used in the context of linear regression. And of course there’s the possibility that I’ve just mangled something else entirely.
Does anything jump out as dodgy or wrong? Or is this even a good approach at all?
