How can I get the credible interval from the output of a Gaussian Process? All I have is the mean and covariance.

When looking at the example for GP in numpyro, I see this line:

`sigma_noise = jnp.sqrt(jnp.clip(jnp.diag(K), a_min=0.0)) * jax.random.normal(rng_key, X_test.shape[:1])`

It looks like it would be possible to make samples by taking the square root of the diagonal in the covariance matrix. And then add random gaussian noise a bunch of times.