Modeling Gaussian process from multiple function observations?

I have ~ 1000 timeseries, each 100 timesteps long.

I’d like to use each timeseries as a function observation of the same Gaussian process. Then, given a new timeseries of just the first 10 timesteps, generate samples from the gaussian process conditioned on the new data.

Is this feasible in pymc3? Are there any examples to help get started?

It looks like the below discourse topic attempted something similar, but the link to jupyter notebook containing an example is now missing and, since that topic, a new MatrixNormal class has been created.

Yes this is possible, there are some codes in the post you mentioned that should give you a good start.

Also, since your observed is generated from the same GP, you dont need MatrixNormal class for it.