Presumably your problem is that your model looks something like
y ~ some_distribution(E[f(theta)], phi);
or even worse
y ~ some_distribution(f(theta), phi);
y_ave = mean(y);
// only y_ave is stored
The former, where your likelihood explicitly depends on an expectation, cannot be implemented directly in Stan unless you can calculate the integral. If the integral is one-dimensional then you can use the
integrate_ode function to compute the integral, or you can wait for a future version of Stan will which have a a tanh-sinh quadrature function. Alternatively you can try to Monte Carlo it as you suggest but it will be extremely slow and it is unlikely that you will get a sufficiently accurate answer (and the gradient will be even worse).
The latter is more tricky because it’s essentially a missing data problem – keeping only a summary statistics loses a huge amount of information. The best approach would be to analytically marginalize over the likelihood to get an exact likelihood for the summary statistic, although this will be hard to do outside of the Gaussian case. Another option is to treat all
N data points as new latent parameters and then approximate a likelihood for
p(y_summary | y_1, ..., y_N). See the “Measurement Model” section of the manual for some examples.