Fitting Gaussian process to set of distributions

So I am trying to use Gaussian processes to model human function learning (in a reinforcement learning -ish setting). Humans are trying to guess the value of some stimulus based on the feedback they have received before for some other stimulus.

i.e. there is $\vec{x} = (x_1,x_2,x_3,x_4)$ observed, a prediction $\tilde{y}$ is made and then they are returned the true value $y = f(\vec{x})$ where f is unknown. The idea being they improve their prediction capabilities over time by implicitly learning the function.

Now, currently we are modelling two different methods for learning by using two different Gaussian Process kernels.

We're fitting to the data using these two models by making a model-based $n+1$th prediction using the first $n$ data points (n sets of $x_i$ and $y$) to predict the participants' response $\tilde{y}$. The log likelihood $\log p(\tilde{y}^{n+1} | \vec{x}^{1:n},y^{1:n},x^n)$ for each of these predictions was then summed and maximised to obtain the fitted model parameters. These models are then tested on a hold-out set to compare which fits best.

Now, the thing I haven't yet mentioned is that as well as a prediction, we ask participants for a confidence estimate. For the sake of simplicity, let's say they give an estimate of their confidence half-range, equal to $k\sigma$ where $k$ varies for each participant.

So, instead of what we do now, we would like to try and obtain better model fitting by fitting to not only their predictions but their confidence estimates too. We're a bit lost on where to start with this and couldn't find anything in the available literature.

If $k$ is known, one idea I had was to simulate dummy values, such that for each $\tilde{y}$ I generate ~100 predictions drawn from the distribution $\mathcal{N}(\tilde{y},\sigma^2)$, and then fit the GP to the whole training set of data using ML. One problem with this is we're then just fitting the best model to the participants predictions, rather than taking into account the learning over time (that we vaguely tackled in our one-step-ahead method).

Other than that, we haven't been able to find any reasonable method to fitting to these confidence estimates, and I thought I would turn to this community for help. Does anyone know of any way of fitting a gaussian process to a set of means AND variances/SDs? Is there such way that might make sense in making 'one-step-ahead' predictions? Any ideas for known or unknown k would be hugely appreciated, even if you're just half-guessing something that might work (I will simulate data using GPs first to see whether any of these methods retrieves the parameters well).

If anything is unclear, please let me know. I understand this post is actually quite light on mathematical content; I guess the main reason being I'm more questioning a concept or idea without being stuck on a particular method itself. I hope this kind of question is okay.

Best,

H

EDIT: After trying to potentially use the sigma calculated from participant responses of confidence as the $\sigma_n$ in the Gaussian Process, unsuccessfully, I'm still looking for a solution to this. Thanks all