TransformedTargetRegressor is incompatible with GaussianProcessRegressor when using return_std=True

### Describe the bug

`TransformedTargetRegressor` is incompatible with `GaussianProcessRegressor` when using `return_std=True` due to the fact that the model returns a tuple rather than an array.



### Steps/Code to Reproduce

```python
import numpy as np

from sklearn.preprocessing import PowerTransformer
from sklearn.compose import TransformedTargetRegressor
from sklearn.gaussian_process import GaussianProcessRegressor

from sklearn.datasets import fetch_openml
from sklearn.utils import shuffle

from sklearn.model_selection import train_test_split
from sklearn.compose import make_column_transformer
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import OrdinalEncoder
from sklearn.compose import make_column_selector

def load_ames_housing():
    df = fetch_openml(name="house_prices", as_frame=True)
    X = df.data
    y = df.target

    features = [
        "YrSold",
        "HeatingQC",
        "Street",
        "YearRemodAdd",
        "Heating",
        "MasVnrType",
        "BsmtUnfSF",
        "Foundation",
        "MasVnrArea",
        "MSSubClass",
        "ExterQual",
        "Condition2",
        "GarageCars",
        "GarageType",
        "OverallQual",
        "TotalBsmtSF",
        "BsmtFinSF1",
        "HouseStyle",
        "MiscFeature",
        "MoSold",
    ]

    X = X[features]
    X, y = shuffle(X, y, random_state=0)

    X = X[:600]
    y = y[:600]
    return X, np.log(y)

X, y = load_ames_housing()

cat_selector = make_column_selector(dtype_include=object)
num_selector = make_column_selector(dtype_include=np.number)

cat_tree_processor = OrdinalEncoder(
    handle_unknown="use_encoded_value", unknown_value=-1
)
num_tree_processor = SimpleImputer(strategy="mean", add_indicator=True)

tree_preprocessor = make_column_transformer(
    (num_tree_processor, num_selector), (cat_tree_processor, cat_selector)
)

X_train, X_test, y_train, y_test = train_test_split(X, y)

X_train = tree_preprocessor.fit_transform(X_train)
X_test = tree_preprocessor.transform(X_test)

gpr = GaussianProcessRegressor(random_state=0)
yj = PowerTransformer("yeo-johnson")
tt_gpr = TransformedTargetRegressor(gpr, transformer=yj)

tt_gpr.fit(X_train, y_train)

tt_gpr.predict(X_test, return_std=True)
```

### Expected Results

Predictions are returned from the transformed regressor. It is potentially unclear how the transformer should apply to the predictions given that custom behavior is likely required for each type of transformer breaking the `func`, `inverse_func` API. However if this is the case it would be good to have a more informative error message/something in the documentation.

### Actual Results

```
AttributeError                            Traceback (most recent call last)
[<ipython-input-3-43ceeb18ba48>](https://localhost:8080/#) in <module>
     31 tt_gpr.fit(X_train, y_train)
     32 
---> 33 tt_gpr.predict(X_test, return_std=True)

[/usr/local/lib/python3.7/dist-packages/sklearn/compose/_target.py](https://localhost:8080/#) in predict(self, X, **predict_params)
    273         check_is_fitted(self)
    274         pred = self.regressor_.predict(X, **predict_params)
--> 275         if pred.ndim == 1:
    276             pred_trans = self.transformer_.inverse_transform(pred.reshape(-1, 1))
    277         else:

AttributeError: 'tuple' object has no attribute 'ndim'
```

### Versions

```shell
System:
    python: 3.7.14 (default, Sep  8 2022, 00:06:44)  [GCC 7.5.0]
executable: /usr/bin/python3
   machine: Linux-5.10.133+-x86_64-with-Ubuntu-18.04-bionic

Python dependencies:
          pip: 21.1.3
   setuptools: 57.4.0
      sklearn: 1.0.2
        numpy: 1.21.6
        scipy: 1.7.3
       Cython: 0.29.32
       pandas: 1.3.5
   matplotlib: 3.2.2
       joblib: 1.2.0
threadpoolctl: 3.1.0

Built with OpenMP: True
```
```

COMMENT:
For anyone with the same issue, wanting to use a TransformedTargetRegressor with predictions that include standard deviations, I have used a wrapper in a project I am working on. Unfortunately, it is not as straightforward as the usual func/inverse_func pattern as we should adjust the mean prediction based on the variance and adjust the standard deviations when inverting to get back to original data scale.

Here is the wrapper:
```

from sklearn.base import BaseEstimator
from sklearn.compose import TransformedTargetRegressor

from ml_app.modelling.entities.numeric_transformers import Transformer
from ml_app.utils import func_supports_arg

class TransformedTargetRegressorWithUncertainty(TransformedTargetRegressor):
    """
        Thin wrapper over sklearn.compose.TransformedTargetRegressor.
        Allows for predict methods that return a tuple of (mean, std) instead of just mean.
        Also allows for corrected inverse_transform methods that take into account the
        target distribution transformation.
    """
    def __init__(self, regressor: BaseEstimator, transformer: Transformer):
        if transformer.inverse_mean_func is None or transformer.inverse_std_func is None:
            raise ValueError(
                "To support predictions with a standard deviation a transformer"
                "must have inverse_mean_func and inverse_std_func functions.")
        self.inverse_mean_func = transformer.inverse_mean_func
        self.inverse_std_func = transformer.inverse_std_func
        super().__init__(regressor, func=transformer.func, 
                         inverse_func=transformer.inverse_func)

    def predict(self, X, return_std=False, with_obs_noise=False, **predict_params):
        """
        Predict using the underlying regressor and transform the result back.
        """
        if return_std:
            model: BaseEstimator = self.regressor_
            if func_supports_arg(model.predict, 'return_std'):
                # this is a parameter I use for a GPy SparseGP wrapper, you can remove if not needed
                if 'with_obs_noise' not in predict_params:
                    predict_params['with_obs_noise'] = with_obs_noise
                tran_pred, tran_std = model.predict(X, return_std=True, **predict_params)
            else:
                # may have to support other scenarios in the future e.g. if they have a 
                # separate predict_std method
                tran_pred, tran_std = model.predict(X, **predict_params)
            pred = self.inverse_mean_func(tran_pred, tran_std)
            std = self.inverse_std_func(tran_pred, tran_std)
            return pred, std
        return super().predict(X, **predict_params)  
```

The numeric transformers referred to here are defined as follows. We have transformation for log and square root, but you could add others as needed:

```
from dataclasses import dataclass
from enum import Enum
from typing import Callable

import numpy as np
import pandas as pd
from scipy.stats import skew


def log_plus_tiny(x):
    """We add a tiny value to allow zero inputs"""
    return np.log(x + 1e-7)


def exp_minus_tiny(x):
    """To invert log_plus_tiny"""
    return np.exp(x) - 1e-7


@dataclass
class Transformer:
    description: str
    func: Callable[[np.ndarray], np.ndarray]
    inverse_func: Callable[[np.ndarray], np.ndarray]
    allows_negatives: bool 
    # these are only used if the model provides standard deviations with predictions
    inverse_mean_func: None|Callable[[np.ndarray, np.ndarray], np.ndarray] = None
    inverse_std_func: None|Callable[[np.ndarray, np.ndarray], np.ndarray] = None


class LogTransformer(Transformer):
    def __init__(self):
        super().__init__(
            description="log(x + 1e-7)",
            func=log_plus_tiny, 
            inverse_func=exp_minus_tiny, 
            allows_negatives=False,
            inverse_mean_func=inverse_log_mean,
            inverse_std_func=inverse_log_std
        )


class SqrtTransformer(Transformer):
    def __init__(self):
        super().__init__(
            description="sqrt(x)",
            func=np.sqrt, 
            inverse_func=np.square, 
            allows_negatives=False,
            inverse_mean_func=inverse_sqrt_mean,
            inverse_std_func=inverse_sqrt_std
        )


class IdentityTransformer(Transformer):
    def __init__(self):
        super().__init__(
            description="",
            func=lambda x: x, 
            inverse_func=lambda x: x, 
            allows_negatives=True
        )


# ⚠️ If you change or remove value add a mapping in converters.py
# Values defined here are used in the UI & API dropdowns
class NumericTransformerOption(Enum):
    na = "No transformation"
    log = "Log"
    sqrt = "Square root"

    def get_transformer(self) -> Transformer:
        if self == NumericTransformerOption.log:
            return LogTransformer()
        elif self == NumericTransformerOption.sqrt:
            return SqrtTransformer()
        else:
            return IdentityTransformer()
        

def get_transformer_with_least_skew(
        y: list[float|int|None] | pd.Series | np.ndarray) -> NumericTransformerOption:
    """
    Apply each available transformation to the data and return the one with the 
    lowest abs skew. 
    """
    try:
        y_array = np.array(y, dtype=float)
    except ValueError:
        return NumericTransformerOption.na

    # Filter out NaN values
    y_array = y_array[~np.isnan(y_array)]

    # Handle the case where all values are NaN or the array is empty
    if y_array.size == 0:
        return NumericTransformerOption.na
    
    any_negative = any(y_array < 0)
    min_skew = abs(skew(y_array))
    best_transformation = NumericTransformerOption.na

    # try each enum value that has a transformation function
    for t in NumericTransformerOption:
        transformer = t.get_transformer()
        if isinstance(transformer, IdentityTransformer):
            continue

        if any_negative and not transformer.allows_negatives:
            continue

        transformed_y = transformer.func(y_array)
        transformed_skew = abs(skew(transformed_y))

        if transformed_skew < min_skew:
            min_skew = transformed_skew
            best_transformation = t

    return best_transformation


# needed for Target transformer when we have std
def inverse_log_mean(log_mean, log_std):
    """
    Calculate the mean in the original scale for a log-normal distribution.

    This function transforms the mean from the logarithmic scale back to the 
    original scale. It corrects for the skewness of the log-normal distribution 
    by incorporating the variance.
    """
    log_var = log_std ** 2    
    return np.exp(log_mean + log_var / 2) - 1e-7


def inverse_log_std(log_mean, log_std):
    """
    Calculate the standard deviation in the original scale for a log-normal distribution.

    This function transforms the standard deviation from the logarithmic scale 
    back to the original scale. It leverages the properties of the log-normal 
    distribution, using the exponential transformation to compute the variance 
    and then obtaining the standard deviation.
    """
    log_var = log_std ** 2
    var = np.exp(2 * log_mean + log_var) * (np.exp(log_var)-1)
    return np.sqrt(var)


def inverse_sqrt_mean(sqrt_mean, sqrt_std):
    """
    This function calculates the mean in the original space. 

    The mean of squared values isn't just the square of the mean of the values 
    because squaring is a non-linear operation. A reasonable approximation is 
    to use the second moment about the mean (i.e., the variance plus the square 
    of the mean) of the sqrt-transformed values.
    """ 
    return sqrt_mean**2 + sqrt_std**2


def inverse_sqrt_std(sqrt_mean, sqrt_std):
    """
    This function calculates the standard deviation in the original space. 

    We need to consider how variance (the square of the standard deviation) 
    transforms. The variance of the squared values can be approximated as the 
    fourth moment about the mean minus the square of the second moment about 
    the mean.
    """
    return np.sqrt(4 * sqrt_mean**2 * sqrt_std**2 + 2 * sqrt_std**4)
```

We test the transformation functions above with:

```
class TestInverseFuncsWithUncertinty(TestCase):

    def test_log_inverse_with_simulation(self):
        # Generating a normal distribution - to simulate original log-normal data that
        # has been log transformed ready for training
        np.random.seed(1)  # for reproducibility
        log_values = np.random.normal(loc=1.4, scale=0.7, size=200000) # large needed for 2dp precision
        log_mean = np.mean(log_values)
        log_std = np.std(log_values)

        values = np.exp(log_values)
        expected_mean = np.mean(values)
        expected_std = np.std(values)

        mean = inverse_log_mean(log_mean, log_std)
        self.assertAlmostEqual(expected_mean, mean, places=2)

        std = inverse_log_std(log_mean, log_std)
        self.assertAlmostEqual(expected_std, std, places=2)

    def test_sqrt_inverse_with_simulation(self):
        # Generating a normal distribution - to simulate original skewed data that
        # has been sqrt transformed ready for training
        np.random.seed(1)  # for reproducibility
        sqrt_values = np.random.normal(loc=1.4, scale=0.7, size=200000) # large needed for 2dp precision
        sqrt_mean = np.mean(sqrt_values)
        sqrt_std = np.std(sqrt_values)

        values = np.square(sqrt_values)
        expected_mean = np.mean(values)
        expected_std = np.std(values)

        mean = inverse_sqrt_mean(sqrt_mean, sqrt_std)
        self.assertAlmostEqual(expected_mean, mean, places=2)

        std = inverse_sqrt_std(sqrt_mean, sqrt_std)
        self.assertAlmostEqual(expected_std, std, places=2)
```

Labels: Enhancement, module:compose