Thank you for your response. We really appreciate it.
* **"...I still think that the paper would be much stronger if the authors demonstrated that the method scales to parameter spaces with dimensionality $>>$ 2":** Agreed. Having an example with more than 4 parameters (as is the case with the Turin model in Section 5) can further strengthen our paper. To that end, we ran our method on the 10-dimensional Gaussian linear example with fixed covariance matrix $\Sigma$ (parameter of interest is the mean vector) used in the RNPE paper and also available in the SBI benchmark library. To make the model misspecified, we used the same contamination model for the observed data used in the paper, i.e., $\mathbb{Q} = (1-\epsilon)\mathcal{N}(\theta\_{\mathrm{true}}, \Sigma) + \epsilon \mathcal{N}(\theta\_c, \Sigma)$ where $\epsilon = 10$%, $\theta\_{\mathrm{true}} = [0.5, \dots, 0.5]^\top$, $\theta\_c = [2,\dots,2]^\top$, $p(\theta) = \mathcal{U}([-1,1])^{10}$. The average MMD between the posterior predictive distribution and the observed data over 100 runs is shown in the following table (std. deviation is reported in the parenthesis). We will of course include these results in the paper.
| | NPE | NPE-RS | NPE-RS | NPE-RS |
|:---:|:----:|:--------------:|:------------:|:-------------:|
| | - | $\lambda = 20$ | $\lambda=50$ | $\lambda=100$ |
| **MMD** | 0.26 (0.02) | 0.19 (0.04) | **0.18** (0.06) | 0.21 (0.08) |
* **"In my opinion, claims like 'It provides accurate inference results even when the model is well-specified' (L 67) are not warranted and have to be removed or significantly trimmed down":** Apologies if we weren't clear before. We agree that our results do not support the claim that "It provides accurate inference results even when the model is well-specified, thus circumventing any need to detect model misspecification", and we will remove it.
* **"Wouldn't negative log-likelihood on misspecified x be a better measure?":** The log-likelihood is not applicable for simulator-based models due to the intractability of the likelihood function. We did consider fitting a multivariate Gaussian distribution to the predictive distributions from Ricker, OUP, and Turin model, and evaluating the log-likelihood of the data under the respective Gaussians. However, that introduces a new level of misspecification, as none of these models produce data that is jointly Gaussian. Moreover, when the observed data is corrupted by a few outlier points (thus causing model misspecification), the sum of log-likelihoods gets dominated by the few outliers, even if most of the observed data is explained well by the model. This is the reason why people have proposed generalised Bayesian inference (GBI) frameworks, where the likelihood term is replaced by a robust loss such as the MMD (which is known to be robust to outliers) to account for misspecification (see Knoblauch et al. (2019) and ref [24] for more details). We therefore used the MMD between the posterior predictive distributions and the observed data (bottom row of Figure 2) as the measure for robustness to misspecification (apart from RMSE).
Knoblauch, J., Jewson, J., \& Damoulas, T. (2019). Generalized variational inference: Three arguments for deriving new posteriors. arXiv preprint arXiv:1904.02063.
* **"Which of the three metrics shown in figure 5 should correspond exactly to "robustness to misspecification" and why is this warranted?":** In Figure 5, we are testing the claim that our method indeed converges to the NPE posterior as $\lambda$ goes to zero, and converges to the prior as $\lambda$ goes to infinity. To do that, we varied $\lambda$ and computed the distance (in terms of MMD) between the NPE-RS posterior and the prior (Fig 5(left)), and between the NPE-RS posterior and the NPE posterior (middle and right).
* **Does this mean that using cross-validation requires several datasets of observed data (which are potentially misspecified and can not easily be synthetically generated) in order to choose $\lambda$? What if only one dataset is available at training time?:** We do not assume that multiple observed datasets are available, even though that is the case in some fields like Astrophysics where available data is plenty. We meant that for one observed dataset $\mathbf{y}\_{1:n} = \{\mathbf{y}^{(1)}\, \dots, \mathbf{y}^{(n)}\}$ with $n$ iid samples available at training time, we can use a random subset of the data, say $\mathbf{y}\_{1:m}$, where $m<n$, as the validation set to choose $\lambda$, and the rest of the $n-m$ points to run NPE-RS (which is what we did for the radio propagation experiment). In the limiting case that only one observed data point is available, there is not going to be any significant posterior update anyway.