## EDIT: I’m dumb and used the wrong curve for reference. Stan actually did great.

For reference:

Thanks for the help so far - I suppose I will be back later as I try to extend this further with the rest of my data.

Hello all,

I followed the approach suggested by @bbbales2. I actually found a better parametrization for the underlying components which requires a single parameter for each one and after some exploration, it seems that only three of the five hypothesized components are present here. Regardless, I’ve been getting some odd results.

The fit that HMC settled for looks quite off:

(the intensity is normalized here and it isn’t the same curve I initially showed, but this one should actually follow the model more closely given how the experiments were conducted)

The results were consistent over a few sets of reasonable initial values. I decided to see if R’s `optim()`

could do better, and it seemingly did:

However, this fit was very susceptible to the choice of initial values.

I thought that using `rstan::optimizing()`

might achieve a compromise of stabiltiy and good fits, but it gave essentially identical results to those from HMC.

Since Stan’s solutions seem obviously suboptimal, I’m wondering if there might be something wrong with my specification. But the unstable fits from `optim()`

makes me think that maybe I need to add some additional restrictions? Not sure how to proceed from here. I’ve got ~60 more curves where the one I’m trying to recover here is present mixed with more components, but if I can’t reliably recover the “pure” signal I don’t suppose using the entire dataset makes sense.

Here’s the model I’ve currently set up (the underlying components are parametrized solely in terms of their maxima):

```
data {
int<lower=0> n; // n# observations
int<lower=0> k; // n# components (tried 2-5, but 3 works best)
vector[k] lm0; // prior maxima
vector[k] s0; // how much to trust priors
vector[n] I; // intensity (measured signal)
vector[n] v; // wavenumber (reciprocal of wavelength shown in graph)
}
parameters {
// Maxima for the components
ordered[k] lm;
// Relative intensities
real<lower=0,upper=1> pu[k];
// Noise
real<lower=0> error;
}
model {
// Intermediate values to ease writing the component formula
vector[k] a;
vector[k] r;
vector[k] vm;
vector[k] vneg;
vector[k] vpos;
// Signals
vector[k] Im[n]; // Latent singal contributed per component
vector[n] mI; // Expected total signal
// Component formula
// Place prior on maxima and invert it
lm ~ normal(lm0,s0);
for(i in 1:k){
vm[i] = (10^7)/lm[i];
}
// Intermediate calculations
vneg = 1.177*vm - 7780;
vpos = 0.831*vm + 7070;
// Per-component emission formula
for(i in 1:n){
mI[i] = 0;
}
for(j in 1:n){
for(i in 1:k){
r[i] = (vm[i] - vneg[i])/(vpos[i] - vm[i]);
a[i] = vm[i] +r[i]*(vpos[i] - vneg[i])/(r[i]^2 - 1);
Im[j][i] = exp(-log(2)*(log(a[i] - v[j])-log(a[i] - vm[i]))^2/(log(r[i])^2));
mI[j] = mI[j] + pu[i]*Im[j][i];
}
}
// Prior for noise
error ~ normal(0, 10);
// Sampling
I ~ normal(mI, error);
}
```