Hi @essicolo, thank you for sharing.
This is a very nice example of the intrinsic model of coregionalization (ICM) using Pymc.
In the model, I am not sure if * is a Kronecker product (⊗) ?
## Combined covariance
cov_f = coreg * cov_feature # (B * K) or (B⊗K)?
There is a pymc.math.kron_dot function in pymc, but I am not sure it will work for this case.