In GaussianProcessRegressor, sampling in the prior before calling fit via sample_y will assume that y is made of a single target. However, this is not necessarily the case. Therefore, the shape of the output of sample_y before and after fit is different.

In order to solve this inconsistency, we need to introduce a new parameter n_targets=None. Before calling fit this parameter should be explicitly set by the user. After fit, we can use the information of the target seen during fit without explicitly setting the parameter.