Support and test in-place sampling

[devmotion]

Currently something like

julia> gp = GP(SqExponentialKernel());

julia> x = rand(11);

julia> rand!(f(x), similar(y))

throws an error since Distributions._rand! is not defined for FiniteGPs (this is the method that Distributions falls back to when one calls rand! for multivariate distributions and that one is supposed to define for new multivariate distributions: https://juliastats.org/Distributions.jl/latest/extends/#Multivariate-Sampler).

This PR adds a definition of Distributions._rand! for FiniteGPs. Distributions also falls back to it when calling rand but I think often this is a bad idea (one has to use type heuristics that fail sometimes - e.g. tests fail if one removes the custom definitions in AbstractGPs) and so I did not remove or modify the existing definitions of rand. It might still be useful to define Distributions._rand! instead of Random.rand! since then e.g. Distributions already checks if the size of the input is correct and one gets support for arrays of vectors for free (see https://github.com/JuliaStats/Distributions.jl/blob/a0cb0969d755872586ec79e985503f03beae9d08/src/multivariates.jl#L47-L66).

[codecov-commenter]

Codecov Report

Merging #176 (1a105a6) into master (0e223c8) will increase coverage by 0.01%.
The diff coverage is 100.00%.

@@            Coverage Diff             @@
##           master     #176      +/-   ##
==========================================
+ Coverage   99.23%   99.24%   +0.01%     
==========================================
  Files           9        9              
  Lines         260      266       +6     
==========================================
+ Hits          258      264       +6     
  Misses          2        2              

Impacted Files	Coverage Δ
src/abstract_gp/finite_gp.jl	100.00% <100.00%> (ø)

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 0e223c8...1a105a6. Read the comment docs.

[theogf]

Should we actually replace this by Distributions._rand! instead?

[devmotion]

I guess this would be reasonable if we want to follow the conventions in Distributions. I am just not sure if we want to do this, it feels weird if an internal method such as _rand! is part of the API (of course, this is a problem of Distributions but since we already define rand separately maybe we want to not follow it here either).

[theogf]

if we want to follow the conventions in Distributions

That seems like something we want to do no? Maybe we can add a reference to their API docs for better understanding?

[devmotion]

But if we really want to follow their conventions we would/should not have to implement rand. However, I really think we should not do this since it is quite likely to lead to incorrect types (I checked and tests fail since the output type is incorrect for Float32 - of course, one could work around this by defining eltype(::FiniteGP) but I think the easier solution is to only perform in-place sampling if requested 🤷). So we do not follow the abstractions in Distributions even if we tell users to implement Distributions._rand!.

[willtebbutt]

The other reason to like having an explicit implementation for rand is AD, since our AD tools really don't like mutation.

[willtebbutt]

I'm also happy with this.