Implement SArray-Compatible Filtering

[THargreaves]

This PR makes some small changes to the filtering interface to allow for full compatibility with SArrays for both closed-form filtering as well as simulation/sampling.

The initial plan we discussed for this approach was to implement a custom Gaussian container. I started this but realised I was running into a lot of edge-cases. Stuff like the computation of the auxiliary statistic (which uses mean on the Gaussian predictive distribution), which would have required us to implement a fair amount of the Distributions.jl interface.

Instead, I thought as a near-term solution, it may be best to work around the limitations of Distributions.jl by making a few changes, rather than trying to reinvent the wheel. In the future, if we decide there are large deviations from the interface that we need and have a clear idea how to implement our own interface, we can pick up that work then.

The primary change of this PR is reverting to using MvNormal for the predictive and filtered distributions of the Kalman filter.

To make this compatible with StaticArrays, we added a layer of abstraction SSMProblems.simulate_from_dist which we can intercept and generate a SArray sample.

The only real consideration with MvNormal is that it requires the covariance to be a PDMat. I.e. it will automatically compute the Cholesky for us which could be inefficient.

It turns out though that this is actually an efficiency gain. I did a FLOPS analysis, and for all dimensions, when computing AΣA', it is faster to first compute the Cholesky of Σ then do B = AL, BB', using PDMats.X_A_Xt than to do the direct computation.

The only issue I'm running into is during the backward recursion of the RTS smoother. Σ_pred - Σ_back is always PD but can be singular if the system is not fully observed (e.g. only observe the position dimension). We don't need to invert this matrix so we could still do something like X_A_Xt but we need to be careful that round-off means the diagonals could become a tiny negative number. I'm not sure of the best way to handle this—just clip to zero?

Any thoughts would be appreciated. I'm not dead set on this approach, so if I custom Gaussian is still preferred, I'm happy to go back to that.

[github-actions]

[JuliaFormatter] _{reported by reviewdog 🐶}

Suggested change

	ys = [SVector{1, T}(y) for y in ys]
	ys = [SVector{1,T}(y) for y in ys]

[github-actions]

SSMProblems.jl/SSMProblems documentation for PR #121 is available at:
https://TuringLang.github.io/SSMProblems.jl/SSMProblems/previews/PR121/