Python / R package for multivariate gaussian process regression for circular data? I have two circular predictor features (two angles between 0 and 360 degrees) and a circular outcome (another angle, between 0 and 360 degrees). I'd like to be able to fit a model and get predictions of the outcome along with their certainty / confidence.
Normally for non-circular data I'd try something like multivariate gaussian process regression (like this tutorial for example). However, I haven't been able to find an appropriate equivalent for circular data. Therefore, I'm looking for any suggestions for how the non-circular case might be done? (in Python / R / Matlab?)
I've seen some suggestions to add/replace(?) the circular predictor features with their cos and sin transformations before running a linear model; although I'm unsure whether this also needs to be done for the circular outcome variable and how that works if its meant to only be one value?

You don't have to transform your outcome variable to model the two non-outcome circular variables, however since your outcome variable is bounded you might have issues with predictions outside the interval [0, 360], which needs to be addressed with a different linear regression variant.
Here is an example in R with three circular variables using a linear model with a sine and cosine term for the two non-outcome variables, showing that you do not need to transform your outcome variable.
set.seed(123)

df=data.frame(
  y=c(sort(runif(50,0,180)),sort(runif(50,180,360))),
  x1=c(sort(runif(50,0,180)),sort(runif(50,180,360))),
  x2=c(sort(runif(50,180,360)),sort(runif(50,0,180)))
)

mod=lm(y~sin(2*pi*x1/360)+cos(2*pi*x1/360)+sin(2*pi*x2/360)+cos(2*pi*x2/360),data=df)

plot(predict(mod))


To get a prediction interval:
predict(mod,interval="predict",level=0.97)
          fit         lwr      upr
1   157.45307  39.6371153 275.2690
2   152.85649  35.2661083 270.4469
3   151.01962  33.2764518 268.7628
4   108.11836 -10.6712443 226.9080
5   104.66850 -13.1176100 222.4546
...

