code: 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):
          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