Here is a really simple model. I’ll just blabber a little. Read what you want ignore what you don’t want.

The idea (people been doing this since at least the 70s) is measuring the natural resonance modes of a lump of metal and trying to estimate the elastic constants.

I get really tight confidence intervals for N as low as 20. What makes me nervous about the inference we’re doing here is that y is really just a list of resonance modes from a single sample of metal. So it’s like we’re doing statistics with a single sample. But! The results we have are very consistent with the reference mechanical data we have (literature) and the measurements we make (Xray diffraction stuff), so I like it (and it works a loooot better than trying to do an optimization).

Now that I type this – I have an unresolved suspicion that if I didn’t model the crystal-sample misorientation (q) is the thing that’s causing this to be hard to estimate. That’s the unit_vector parameter. There’s a ton of symmetry in the system, so technically the posterior on that should be super-multimodal, but usually it isn’t too radical (maybe two modes max?). I need to test if this causes the sampling difficulties instead of just talking about it…

The Xtal-sample misorientation parameter (q) is the misorientation (as a quaternion) between the sample reference frame and the Xtal lattice. This is important for the resonance modes. Think of rocksalt, where your atoms are arranged in a cubic structure. The misorientation parameter means that cubic grid of atoms (xtal reference frame) doesn’t need to line up with the actual cubes of salt (the sample reference frame).

One of the big issues with this is I really don’t want to turn down the precision of the resonance mode calculations (lower precision calculations is one way I can do my inferences faster). I know this can make weird biases if the model is off (I’ve seen it). I think we got really lucky with our data/model in that it actually works really well, but I gotta keep P high (order of Rayleigh-Ritz approximation to eigenvalue problem) or I might accidentally be biasing my results (and Stan won’t be able to help me there).

```
// These are some externally defined functions that do the RUS forward model.
// That is, they take in elastic constants and produce resonance modes
functions {
vector mech_init(int P, real X, real Y, real Z, real density);
matrix mech_rotate(matrix C, vector q);
vector mech_rus(int P, int N, vector lookup, matrix C);
}
// Input data
data {
int<lower = 1> P; // Order of polynomials for Rayleigh-Ritz approx
int<lower = 1> L; // This is a function of P :(
int<lower = 1> N; // Number of resonance modes
// Sample known quantities
real<lower = 0.0> X;
real<lower = 0.0> Y;
real<lower = 0.0> Z;
real<lower = 0.0> density;
// Resonance modes
vector[N] y;
}
transformed data {
vector[L * L * 3 * 3 + L * L + L * L * 3 * 3 * 21 + L * L * 3 * 3] lookup;
lookup = mech_init(P, X, Y, Z, density);
}
// Parameters to estimate
parameters {
real<lower = 0.0> c11;
real<lower = 0.0> a;
real<lower = 0.0> c44;
real<lower = 0.0> sigma;
unit_vector[4] q; // rotation between sample & xtal axes
}
// Build a 6x6 stiffness matrix and rotate it
transformed parameters {
real c12;
matrix[6, 6] C;
c12 = -(c44 * 2.0 / a - c11);
for (i in 1:6)
for (j in 1:6)
C[i, j] = 0.0;
C[1, 1] = c11;
C[2, 2] = c11;
C[3, 3] = c11;
C[4, 4] = c44;
C[5, 5] = c44;
C[6, 6] = c44;
C[1, 2] = c12;
C[1, 3] = c12;
C[2, 3] = c12;
C[3, 2] = c12;
C[2, 1] = c12;
C[3, 1] = c12;
C = mech_rotate(C, q);
}
model {
// Specify a prior on noise level. Units are khz, we expect ~100-300hz in a good fit
sigma ~ normal(0, 2.0);
// Specify soft priors on the elastic constants
c11 ~ normal(1.0, 2.0);
a ~ normal(1.0, 2.0);
c44 ~ normal(1.0, 2.0);
// Resonance modes are normally distributed around what you'd expect from an RUS calculation
y ~ normal(mech_rus(P, N, lookup, C), sigma);
}
```