def analytical_linearisation(self, m, hyp=None):
        """
        Compute the Jacobian of the state space observation model w.r.t. the
        function fₙ and the noise term rₙ.
        The implicit observation model is:
            h(fₙ,rₙ) = E[yₙ|fₙ] + √Var[yₙ|fₙ] rₙ
        The Jacobians are evaluated at the means, fₙ=m, rₙ=0, to be used during
        extended Kalman filtering and extended Kalman EP.
        """
        Jf, Jr = jacrev(self.observation_model, argnums=(0, 1))(m, 0.0, hyp)
        return Jf, Jr