C++ sample::cont_params方法代码示例

    sample transition(sample& init_sample) {
        this->sample_stepsize();

        this->seed(init_sample.cont_params());

        this->hamiltonian_.sample_p(this->z_, this->rand_int_);
        this->hamiltonian_.init(this->z_);

        ps_point z_init(this->z_);

        double H0 = this->hamiltonian_.H(this->z_);

        for (int i = 0; i < L_; ++i)
            this->integrator_.evolve(this->z_, this->hamiltonian_,
                                     this->epsilon_);

        double h = this->hamiltonian_.H(this->z_);
        if (boost::math::isnan(h)) h = std::numeric_limits<double>::infinity();

        double acceptProb = std::exp(H0 - h);

        if (acceptProb < 1 && this->rand_uniform_() > acceptProb)
            this->z_.ps_point::operator=(z_init);

        acceptProb = acceptProb > 1 ? 1 : acceptProb;

        return sample(this->z_.q, - this->hamiltonian_.V(this->z_), acceptProb);
    }

      sample
      transition(sample& init_sample,
                 interface_callbacks::writer::base_writer& info_writer,
                 interface_callbacks::writer::base_writer& error_writer) {
        // Initialize the algorithm
        this->sample_stepsize();

        nuts_util util;

        this->seed(init_sample.cont_params());

        this->hamiltonian_.sample_p(this->z_, this->rand_int_);
        this->hamiltonian_.init(this->z_, info_writer, error_writer);

        ps_point z_plus(this->z_);
        ps_point z_minus(z_plus);

        ps_point z_sample(z_plus);
        ps_point z_propose(z_plus);

        int n_cont = init_sample.cont_params().size();

        Eigen::VectorXd rho_init = this->z_.p;
        Eigen::VectorXd rho_plus(n_cont); rho_plus.setZero();
        Eigen::VectorXd rho_minus(n_cont); rho_minus.setZero();

        util.H0 = this->hamiltonian_.H(this->z_);

        // Sample the slice variable
        util.log_u = std::log(this->rand_uniform_());

        // Build a balanced binary tree until the NUTS criterion fails
        util.criterion = true;
        int n_valid = 0;

        this->depth_ = 0;
        this->divergent_ = 0;

        util.n_tree = 0;
        util.sum_prob = 0;

        while (util.criterion && (this->depth_ <= this->max_depth_)) {
          // Randomly sample a direction in time
          ps_point* z = 0;
          Eigen::VectorXd* rho = 0;

          if (this->rand_uniform_() > 0.5) {
            z = &z_plus;
            rho = &rho_plus;
            util.sign = 1;
          } else {
            z = &z_minus;
            rho = &rho_minus;
            util.sign = -1;
          }

          // And build a new subtree in that direction
          this->z_.ps_point::operator=(*z);

          int n_valid_subtree = build_tree(depth_, *rho, 0, z_propose, util,
                                           info_writer, error_writer);
          ++(this->depth_);

          *z = this->z_;

          // Metropolis-Hastings sample the fresh subtree
          if (!util.criterion)
            break;

          double subtree_prob = 0;

          if (n_valid) {
            subtree_prob = static_cast<double>(n_valid_subtree) /
              static_cast<double>(n_valid);
          } else {
            subtree_prob = n_valid_subtree ? 1 : 0;
          }

          if (this->rand_uniform_() < subtree_prob)
            z_sample = z_propose;

          n_valid += n_valid_subtree;

          // Check validity of completed tree
          this->z_.ps_point::operator=(z_plus);
          Eigen::VectorXd delta_rho = rho_minus + rho_init + rho_plus;

          util.criterion = compute_criterion(z_minus, this->z_, delta_rho);
        }

        this->n_leapfrog_ = util.n_tree;

        double accept_prob = util.sum_prob / static_cast<double>(util.n_tree);

        this->z_.ps_point::operator=(z_sample);
        this->energy_ = this->hamiltonian_.H(this->z_);
        return sample(this->z_.q, - this->z_.V, accept_prob);
      }

 sample transition(sample& init_sample) {
   this->seed(init_sample.cont_params(), init_sample.disc_params());
   return sample(this->_z.q, this->_z.r, - this->_hamiltonian.V(this->_z), 0);
 }