Goodness of fit - stochastic model Going back to a very old question... How to quantify the goodness of fit of a model. Now, I now this is a very general questions, that is probably still the subject of research nowadays, so I will present a more specific case and my initial solution.
I have an experiment where I measure two different quantities ($X$ and $Y$) for different timepoints. This two quantities have very different scales. At each timepoint I have several repetitions of said measurements. This repetitions included both technical and biological replicates. As I expect (and observed) the biological variation to be larger than the technical variation, I end up summarising the data with the mean of each quantity per timepoint and the standard error of the mean (SEM) for each (calculated as sd(data)/sqrt(size of data)).
Now, the model part. My model is a non-linear model whit either 3, 4 or 5 parameters (I'm testing different models, but that's not the point). The model is stochastic, so for every parameter set I summarise the model output by running a specific number of simulations (30) and compute the mean and SEM for each variable in each timepoint. 
After I've done with computing the best set of parameters for the model (another yet related problem) I want to assess the quality of the fitting... Any ideas?


*

*The simplest is that I compute some $r^2$, either weighted by the error in the data or not.

*In my optimization process I computed an analogous to the Likelihood of the data given the parameters:
$L(data | parameters) = \int N(x | u_p, \sigma_p)*N(x | u_o, \sigma_o) dx $ with $u_p$ and $\sigma_p$ being the predicted mean and SEM from the model and $u_o$ and $\sigma_o$ being the observed.
 A: After some time and without any other comments or discussion, I decided to post my proposed solution. Bare in mind that I'm not sure if this is the best solution, but I would welcome any comments or corrections on it.
We start by assuming that both the data measurement and the model simulations behave like a normal distribution. Also, after my fitting procedure I obtain a posterior distribution of the parameters (not just a point estimate).
 I sample $n$ parameters sets $\theta_j$ from the posterior distribution and assume that, for every one, I run enough simulations that I'm confidant on the mean value expected from that specific parameter combination $j$.
I compute the sum of standardize error $Z$ as
$Z = \sum_{h=1}^n \sum_{g=1}^q(\dfrac{(p_{gn}-\tilde d_g)}{\sigma d_g})^2$
where  $n$ is the number of sampled parameter sets, $q$ is the number of datapoints, $p_{gn}$ is the predicted value for timepoint $g$ with $\theta_n$, $\tilde d_g$ is the mean observed value at timepoint $g$ and $\sigma d_g$ is the observed variance at timepoint $g$.
Now, I consider $Z$ to be the test statistic for a $\chi^2$ test with $n\times q$ degrees of freedom.
