Thanks for the resources! I think I was able to define a custom mean function and infer theta, I will post my code here soon. However, I would still like to figure out how to solve the problem without the gaussian process library (option 2) and compare the results to option 1. Since they should result in the same answer, it would be nice to confirm this.