I tried to estimate parameters of a coupled oscillator network for predicting dyadic rhythmic entrainment, using ODE-regression in Stan.
The text is planned as a possible contribution to the Stan forum.
But before doing that I would like to invite your comments in view of preparing a final version.
That’s a really cool problem and overall, this presentation looks great. Especially the discussion at the end. It reminds me of the same kind of synchronization we get in spoken dialogue between human agents.
In terms of specifics, I didn’t see the point of the overall (a), (b), (c) description of the process. It seems to combine simulation with modeling and inference. I think (though I’m not sure) that you’re just trying to lay out the standard paradigm of Bayesian inference.
To report on Stan, I suggest reporting which interface you used (CmdStanPy, RStan, etc.).
That’s a big computer on which to run Stan—did you experiment to see how many parallel threads you could run or did you parallelize the likelihood ODEs? That’s where we’ve had nearly linear speedups with number of threads because there’s so much more computation than communication.
In writing out the ODE system, I didn’t understand the dot notation. Were those zeroes? And if so, why isn’t there a dot in the upper left corner?
I’ve never understood violin plots and the liberties they seem to take with data ranges. For example, how can the error be less than zero in figure 6?
I would suggest trying to plot the data in Table 1 or drop it or maybe move it to an appendix.
I would also suggest including the Stan code and scripts you use to run it. That can just be a link to the source code if it’s in a public repo.
In the discussion, there’s a typo: performans → performs. And I’d qualify the comment on a parameter’s variance to be “posterior variance”. It will also have a prior variance.
P.S. “Stan” isn’t an acronym, so we don’t usually capitalize it.
Thanks a lot for you helpful comments.
A new version can be found here:
- I deleted the (a), (b), (c) description
- the interface used is rstan, this is now added to the text
- I parallelize the 40 regression models of a for loop, using the R-package doParallel, using 4 cores per model. In one of the coming weeks I’ll do a measurement of the time it takes for each model, and in global. In earlier tests (with less constrained defined parameters) it happened that a model did not converge while all other models had converged, or that one model took a very long time compared to the other models. The problem is that the server is occupied until the last model is calculated. Otherwise, one has to interrupt the system. However, I am not an expert in parallel computing - there may be better ways to handle this.
- the matrices K and D were badly represented and this has been redone. The matrices create the ODE-system description and so a dot just means: no creation of coupling or phase delay.
- I used box plots instead of violin plots. The errors depicted show the difference of the mean of the posterior distribution of the parameter (fitted) with the defined parameter value (defined). the fitted value can be higher or lower than the defined value. In the end one could take the absolute value to specify the error but I found it useful to report plus and minus in view of a possible bias.
- I moved Table 1 to Appendix B
- I added the Stan code to Appendix A. A didactical script in R showing how to use the code would be useful. I am working on it.
Thanks a lot and if you find ways to improve the Stan code I will be happy to try it.