I think you’re trying to fit what is typically called a non-linear least squares model, not a mixture model.
That should be straightforward in Stan. Something like (untested):
functions {
real phi(real x, real mu, real sigma){
return 1 / (sigma * sqrt(2 * pi())) * exp(-(x - mu)^2/(2 * sigma^2));
}
data {
...
}
parameters {
ordered[2] m;
real<lower=0> s[2];
real<lower=0> sigma;
}
model {
// Priors
...
// Likelihood
for(i in 1:N){
y[i] ~ normal(phi(x[i], m[1], s[1]) + phi(x[i], m[2], s[2]), sigma);
}
}
There’s a lot you can do to speed this up (vectorizing phi() is the obvious thing) and I’d imagine the posterior is a bit difficult, but this should get you started.