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?
 A: 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
...

