docs: def image_position_stochastic(self, source_x, source_y, kwargs_lens, search_window=10,
                                    precision_limit=10**(-10), arrival_time_sort=True, x_center=0,
                                    y_center=0, num_random=1000, verbose=False):
          """
          Solves the lens equation stochastically with the scipy minimization routine on the quadratic distance between
          the backwards ray-shooted proposed image position and the source position.
          Credits to Giulia Pagano
  
          :param source_x: source position
          :param source_y: source position
          :param kwargs_lens: lens model list of keyword arguments
          :param search_window: angular size of search window
          :param precision_limit: limit required on the precision in the source plane
          :param arrival_time_sort: bool, if True sorts according to arrival time
          :param x_center: center of search window
          :param y_center: center of search window
          :param num_random: number of random starting points of the non-linear solver in the search window
          :param verbose: bool, if True, prints performance information
          :return: x_image, y_image
          """
          x_solve, y_solve = [], []
          for i in range(num_random):
              x_init = np.random.uniform(-search_window / 2., search_window / 2) + x_center
              y_init = np.random.uniform(-search_window / 2., search_window / 2) + y_center
              xinitial = np.array([x_init, y_init])
              result = minimize(self._root, xinitial, args=(kwargs_lens, source_x, source_y), tol=precision_limit ** 2, method='Nelder-Mead')
              if self._root(result.x, kwargs_lens, source_x, source_y) < precision_limit**2:
                  x_solve.append(result.x[0])
                  y_solve.append(result.x[1])
  
          x_mins, y_mins = image_util.findOverlap(x_solve, y_solve, precision_limit)
          if arrival_time_sort is True:
              x_mins, y_mins = self.sort_arrival_times(x_mins, y_mins, kwargs_lens)
          return x_mins, y_mins