title: xarray.Dataset.polyfit: fit two variables
body: <p>I'm trying to use xarray.Dataset.polyfit to get the slope &amp; intercept for weight ~ velocity, for each group of exercise and date.</p>
  <p>What seems to happen is that I'm getting the coefficients for set_number and each variable.</p>
  <p>Here's what I'm doing:</p>
  <pre><code>exercise = ['squat', 'deadlift']
  workout_date = pd.date_range('2022-01-01', freq = 'D', periods = 10)
  set_number = np.arange(10)
  
  n_exercises = len(exercise)
  n_workouts = len(workout_date)
  n_sets = len(set_number)
  
  weight = np.array([[20 + set_number*10] * n_workouts] * n_exercises)
  velocity = np.random.randn(n_exercises, n_workouts, n_sets)
  
  ds = xr.Dataset(data_vars = dict(velocity = (['exercise', 'workout_date', 'set_number'], velocity),
                                   weight = (['exercise', 'workout_date', 'set_number'], weight),
                                   ),
                  coords = dict(exercise = exercise,
                                workout_date = workout_date,
                                set_number = set_number
                                )
                  )
  
  ds.polyfit(dim = 'set_number', deg = 1)
  </code></pre>
  <p>I expect to see a new variable for each exercise &amp; workout_date, containing slope &amp; intercept.</p>