>- It could be the example to compare with the true posterior.
To do this, we have now implemented the Ensemble Kalman Smoother with an advection-diffusion transition model matching the data-generating process. We used $10^4$ ensemble members (the maximum feasible on our machine). Fixing the velocity and diffusion parameters to the true values, we obtain the following results:
$MAE_{\mu}$ = 0.0512$\tiny{\pm 0.0019}$, $RMSE_{\mu}$ = 0.0654$\tiny{\pm 0.0028}$, $MAE_{\sigma}$ = 0.0041$\tiny{\pm 0.0000}$, $RMSE_{\sigma}$ = 0.0045$\tiny{\pm 0.0000}$, $CRPS$ = 0.1025$\tiny{\pm 0.0029}$
As expected, the estimated posterior is very close to the true posterior.
For a more appropriate comparison with our approach, we used a state augmentation approach to estimate the velocity and diffusion parameters jointly with the system states. In contrast to the ST-DGMRF approach, we consider initial and transition noise parameters to be fixed in order to avoid divergence of the EnKS. This yields the following results:
$MAE_{\mu}$ = 0.1249$\tiny{\pm 0.1423}$, $RMSE_{\mu}$ = 0.1925$\tiny{\pm 0.2504}$, $MAE_{\sigma}$ = 0.0046$\tiny{\pm 0.0011}$, $RMSE_{\sigma}$ = 0.0061$\tiny{\pm 0.0031}$, $CRPS$ = 0.1624$\tiny{\pm 0.1226}$
While $MAE_{\sigma}, RMSE_{\sigma}$ are lower than with our approach, $MAE_{\mu}, RMSE_{\mu}$ are slightly higher than with the ST-DGMRF model with $L_{\text{temporal}}=4$.
Note that $MAE_{\sigma}$ and $RMSE_{\sigma}$ are expected to increase as well when considering the noise parameters unknown.
Reducing the ensemble size by half, increases $MAE_{\mu}$ and $RMSE_{\mu}$ further to 0.1647$\tiny{\pm} 0.1274$ and 0.2244$\tiny{\pm 0.1765}$ respectively.
>- Missing observations is not a must for every example, is it?
You are right, missing observations are not a must. However, settings with partially observed system states are a focus of our paper and one of the main motivations for developing our approach (see lines 15-25), which is why we designed our experiments accordingly.
>- There are variational methods learning the, e.g. [Frigola 2014](https://papers.nips.cc/paper_files/paper/2014/hash/139f0874f2ded2e41b0393c4ac5644f7-Abstract.html) and [Naesseth 2017](https://arxiv.org/abs/1705.11140). Dual Kalman filter could also learn the dynamical model.
Thank you for sharing these very interesting papers with us. Using Gaussian processes in combination with variational inference, as done in Frigola 2014, is indeed a very promising approach. However, their application seems (so far) limited to small systems (e.g. one variable and its derivative). Applying such an approach to large systems with interacting components would require designing a suitable GP transition with both high-dimensional input and high-dimensional output. This seems beyond the scope of this paper.
>Fig 1. bottom. Three methods perfectly overlap for input data?
Yes, this is expected because the input data points $\mathbf{y}$ are given to all these methods at inference time. Unless the observation noise is set very high, the estimated posterior $p(\mathbf{x}\mid\mathbf{y})$ will thus be very close to the observed values at these points for all considered methods.
If you have any further questions, we would be happy to answer them.