def profile_likelihood(likelihood, *dims):
    """
    Compute the profile (log)likelihood for
    one or more parameters
    
    Parameters
    ----------
    likelihood : n-d array_like
        (log)likelihood
    grid_steps : 1-d array_like
        the resolution (step length) of each
        dimension of the regular grid
    *dims : int, int, ...
        parameters (grid dimensions) for which
        to compute the profile (log)likelihood
    
    Returns
    -------
    1-d or n-d array
        values of the profile (log)likelihood
        on the regular grid. The number of
        dimensions depends on the number
        of `*dims` arguments given.
    """
    Lp = likelihood.copy()
    ax = 0
    
    for d in range(likelihood.ndim):
        if d in dims:
            ax += 1
            continue
        Lp = Lp.max(axis=ax)
    
    return Lp