If you’re comfortable experimenting a little, I’ve written a package that is intended to make this sort of wrangling easier. It’s here: GitHub - JohnGoertz/Gumbi: Gaussian Process Model Building Interface. It uses pymc3 under the hood, but aims to take care of all the reshaping, transforming, and standardization needed.

I’ve been working on it for a while, but I’ve just started making it public, uploading it to PyPI, hosting the documentation on ReadTheDocs, etc. I feel it’s well documented, but you need to compile the documentation locally at the moment. You can also just take a look at the few example notebooks: Gumbi/docs/source/notebooks at main · JohnGoertz/Gumbi · GitHub

In your case, you would store your data as a single tall DataFrame `df`

with one column for “pixel” (0-1007, or whatever is appropriate), one column for “camera” (‘A’ or ‘B’), one column for “channel” (‘R’/‘G’/‘B’), and one for “value” (probably best to normalize between 0 and 1, noninclusive).

The steps with Gumbi would then just be:

```
ds = gmb.DataSet(df, outputs=['value'], logit_vars=['value']) # Omit logit_vars if your data is not normalized to (0, 1)
gp = gmb.GP(ds)
# You may need to pass `sparse=True` for your 6k datapoints. 100 inducing points are the default, otherwise specify with `n_u`.
gp.fit(outputs=['value'], continuous_dims=['pixel'], categorical_dims=['camera', 'channel'], sparse=True)
# Make predictions for each combination of camera/channel and plot
X = gp.prepare_grid()
axs = plt.subplots(2,3, figsize=(18, 8))[1]
for row, camera in zip(axs.T, ['A','B']):
for ax, channel in zip(row, ['R','G','B']):
y = gp.predict_grid(categorical_levels={'camera': camera, 'channel': channel}, with_noise=False)
gmb.ParrayPlotter(X, y).plot(ax=ax)
```

That setup will use a coregionalization kernel to learn correlations between all combinations of camera+channel, inspired by this other implementation by @bwengals.

The package is still in development, so let me know if you need any help or have any suggestions for improvement!