def _find_reasonable_stepsize(self, position, stepsize_app=1):
        """
        Method for choosing initial value of stepsize

        References
        -----------
        Matthew D. Hoffman, Andrew Gelman, The No-U-Turn Sampler: Adaptively
        Setting Path Lengths in Hamiltonian Monte Carlo. Journal of
        Machine Learning Research 15 (2014) 1351-1381
        Algorithm 4 : Heuristic for choosing an initial value of epsilon
        """
        # momentum = N(0, I)
        momentum = np.reshape(np.random.normal(0, 1, len(position)), position.shape)

        # Take a single step in time
        position_bar, momentum_bar, _ =\
            self.simulate_dynamics(self.model, position, momentum,
                                   stepsize_app, self.grad_log_pdf).get_proposed_values()

        acceptance_prob = self._acceptance_prob(position, position_bar, momentum, momentum_bar)

        # a = 2I[acceptance_prob] -1
        a = 2 * (acceptance_prob > 0.5) - 1

        condition = self._get_condition(acceptance_prob, a)

        while condition:
            stepsize_app = (2 ** a) * stepsize_app

            position_bar, momentum_bar, _ =\
                self.simulate_dynamics(self.model, position, momentum,
                                       stepsize_app, self.grad_log_pdf).get_proposed_values()

            acceptance_prob = self._acceptance_prob(position, position_bar, momentum, momentum_bar)

            condition = self._get_condition(acceptance_prob, a)

        return stepsize_app