I very much enjoyed reading the paper by Lin et al. This is an extremely impressive and thorough piece of work. The associated Python framework is extremely well done, and I was very impressed to see so many options given for how the models are fitted and plotted. My only real concern about this work is the functionality for those who are not experts in Python (or Pyro) as the vast majority of users would be. There are very thorough notebooks in the tutorial section of the Github repo but these are not really helpful for those who want to do a quick straightforward model. My guess is that most users want to fit a GP model using the default values for time uncertainty, error variance, kernel choice, etc, and would like a simple guide for how to get their data from the Excel spreadsheet through the PalaeoSTEHM pipeline. Going even further I would strongly encourage the authors to create a proper Python package to simplify the instructions and coding. 

One notable thing I couldn’t see in the notebooks I ran (or in the paper) was convergence checking for the model. I would say this is absolutely vital for having any faith in the results. Models of this complexity can be extremely difficult to obtain convergence on, and there is a whole range of summary stats available for this in Pyro. It should be part of the default workflow everywhere. Perhaps it is for some scripts and I’ve missed it, but it certainly isn’t discussed in the paper. On a similar vein I’d like to know if the model is calibrated, via scoring rules or even just some posterior predictive distributions (though these can be tricky with bivariate uncertainties), and some kind of out-of-sample performance metrics as would be common in standard ML pipelines. 

Otherwise I really enjoyed the paper and I’m super excited to see how this develops. There were a number of poor sentences which I’ve highlighted below but I don’t think I’ve got all of them. It just needs another language check. 

L32: Change-points

L35: I’d put the reference to the GitHub repo here so people can start coding without needing to read the whole paper. 

Table 1: I’d just review some of these definitions. The conditional probability one about conditioning on an unknown quantity doesn’t read quite right. I also think you: should include one for parameter itself; adjust the line spacing or add horizontal lines to separate the entries better; re-write the likelihood one; and change the uncertainty one which seems to be a frequentist definition. 

Fig 1: External

Fig 1 caption: Platforms

L174: I assume the number of change points is fixed and not learnt?

L176: delete ‘and’

L187: Mu(t) or mu(X) (as used in Eq 10)?

L201: will be shown in Section 2

Table 2: I got confused by what the sampling covariance is and how it is calculated for deterministic models. Please expand in the text

Eq 15 (and perhaps others). The usual way to present normal distributions is mean and variance, not sd. 

Fig 5 (bottom right) and others. It always bugs me slightly that the uncertainty in the rate for the present is more unknown when it is the period when we have the most data. Is there a way to solve this with these models? It strikes me that we should be using temporally non-stationary models that allow for far reduced variance (and hence variance on the derivative) to capture the rate of the most recent periods.