<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def approx_fprime(x, f, epsilon=None, args=(), kwargs=None, centered=True):
 """ Gradient of function, or Jacobian if function fun returns 1d array Parameters x : array parameters at which the derivative is evaluated fun : function `fun(*((x,)+args), **kwargs)` returning either one value or 1d array epsilon : float, optional Stepsize, if None, optimal stepsize is used. This is _EPS**(1/2)*x for `centered` == False and _EPS**(1/3)*x for `centered` == True. args : tuple Tuple of additional arguments for function `fun`. kwargs : dict Dictionary of additional keyword arguments for function `fun`. centered : bool Whether central difference should be returned. If not, does forward differencing. Returns ------- grad : array gradient or Jacobian Notes ----- If fun returns a 1d array, it returns a Jacobian. If a 2d array is returned by fun (e.g., with a value for each observation), it returns a 3d array with the Jacobian of each observation with shape xk x nobs x xk. I.e., the Jacobian of the first observation would be [:, 0, :] """

    kwargs = {} if kwargs is None else kwargs
    n = len(x)
    f0 = f(*(x,) + args, **kwargs)
    dim = np.atleast_1d(f0).shape  # it could be a scalar
    grad = np.zeros((n,) + dim, float)
    ei = np.zeros(np.shape(x), float)
    if not centered:
        epsilon = _get_epsilon(x, 2, epsilon, n)
        for k in range(n):
            ei[k] = epsilon[k]
            grad[k, :] = (f(*(x + ei,) + args, **kwargs) - f0) / epsilon[k]
            ei[k] = 0.0
    else:
        epsilon = _get_epsilon(x, 3, epsilon, n) / 2.
        for k in range(n):
            ei[k] = epsilon[k]
            grad[k, :] = (f(*(x + ei,) + args, **kwargs) -
                          f(*(x - ei,) + args, **kwargs)) / (2 * epsilon[k])
            ei[k] = 0.0
    grad = grad.squeeze()
    axes = [0, 1, 2][:grad.ndim]
    axes[:2] = axes[1::-1]
    return np.transpose(grad, axes=axes).squeeze()