If I have access to multiple samples from a Gaussian process with known covariance kernel but unknown parameters (i.e. unknown lengthscale), it is straightforward to estimate the lengthscale using maximum likelihood methods. What I would like to know is how to calculate the full posterior of the lengthscale parameter.
There are various approximations I could use, such as approximate Bayesian computation. But I feel sure there must be a way to do it in closed form.