input: Write a python function to perform multiple linear regression and return the coefficients.
  assert np.allclose(multiple_linear_regression(np.array([[1, 2], [3, 4], [5, 6]]), np.array([1, 2, 3])), [-1.00000000e+00, 5.00000000e-01, 3.62330302e-16])
  assert np.allclose(multiple_linear_regression(np.array([[1, 1], [2, 2], [3, 3]]), np.array([2, 4, 6])), [2.00000000e+00, 4.44089210e-16, 2.00000000e+00])
  assert np.allclose(multiple_linear_regression(np.array([[1], [2], [3], [4]]), np.array([2, 2.5, 3.5, 5])), [1.25 , 0.875])
output: def multiple_linear_regression(X, y):
      import numpy as np
      # Add a column of ones to the input features matrix to account for the intercept
      X = np.column_stack([np.ones(X.shape[0]), X])
      # Calculate the coefficients using the normal equation
      beta = np.linalg.inv(X.T @ X) @ X.T @ y
      return beta