diff --git a/rstan/rstan/NAMESPACE b/rstan/rstan/NAMESPACE index 729054221..3488228d9 100644 --- a/rstan/rstan/NAMESPACE +++ b/rstan/rstan/NAMESPACE @@ -21,7 +21,7 @@ export( stan_model, stanc, stanc_builder, - stan_version, + stan_version, stan, stan_rdump, read_rdump, @@ -35,12 +35,12 @@ export( rstan_options, As.mcmc.list, set_cppo, - stan_plot, - stan_trace, - stan_hist, - stan_dens, + stan_plot, + stan_trace, + stan_hist, + stan_dens, stan_scat, - stan_ac, + stan_ac, stan_diag, stan_par, stan_rhat, @@ -54,8 +54,8 @@ export( cpp_object_initializer, # get_rstan.options check_hmc_diagnostics, - check_divergences, - check_treedepth, + check_divergences, + check_treedepth, check_energy, get_divergent_iterations, get_max_treedepth_iterations, @@ -70,21 +70,21 @@ export( ) exportClasses( - stanmodel, stanfit -) + stanmodel, stanfit +) exportMethods( -# print, plot, -# extract, +# print, plot, +# extract, optimizing, vb, - get_cppcode, get_cxxflags, # for stanmodel + get_cppcode, get_cxxflags, # for stanmodel show, sampling, summary, extract, traceplot, plot, get_stancode, get_inits, get_seed, get_cppo_mode, - log_prob, grad_log_prob, + log_prob, grad_log_prob, hessian_log_prob, hessian_times_vector_log_prob, unconstrain_pars, constrain_pars, get_num_upars, get_seeds, get_adaptation_info, get_sampler_params, - get_logposterior, + get_logposterior, get_posterior_mean, get_elapsed_time, get_stanmodel, diff --git a/rstan/rstan/R/stanfit-class.R b/rstan/rstan/R/stanfit-class.R index 8ef4c4e08..f7af16cfe 100644 --- a/rstan/rstan/R/stanfit-class.R +++ b/rstan/rstan/R/stanfit-class.R @@ -667,6 +667,30 @@ setMethod("grad_log_prob", signature = "stanfit", object@.MISC$stan_fit_instance$grad_log_prob(upars, adjust_transform) }) +if (!isGeneric("hessian_log_prob")) { + setGeneric(name = "hessian_log_prob", + def = function(object, ...) { standardGeneric("hessian_log_prob") }) +} + +setMethod("hessian_log_prob", signature = "stanfit", + function(object, upars, adjust_transform = TRUE) { + if (!is_sfinstance_valid(object)) + stop("the model object is not created or not valid") + object@.MISC$stan_fit_instance$hessian_log_prob(upars, adjust_transform) + }) + +if (!isGeneric("hessian_times_vector_log_prob")) { + setGeneric(name = "hessian_times_vector_log_prob", + def = function(object, ...) { standardGeneric("hessian_times_vector_log_prob") }) +} + +setMethod("hessian_times_vector_log_prob", signature = "stanfit", + function(object, upars, v, adjust_transform = TRUE) { + if (!is_sfinstance_valid(object)) + stop("the model object is not created or not valid") + object@.MISC$stan_fit_instance$hessian_times_vector_log_prob(upars, v, adjust_transform) + }) + setMethod("traceplot", signature = "stanfit", function(object, pars, include = TRUE, unconstrain = FALSE, inc_warmup = FALSE, window = NULL, nrow = NULL, ncol = NULL, diff --git a/rstan/rstan/inst/include/rstan/rcpp_module_def_for_rstan.hpp b/rstan/rstan/inst/include/rstan/rcpp_module_def_for_rstan.hpp index d3f58b996..50e412272 100644 --- a/rstan/rstan/inst/include/rstan/rcpp_module_def_for_rstan.hpp +++ b/rstan/rstan/inst/include/rstan/rcpp_module_def_for_rstan.hpp @@ -25,6 +25,10 @@ RCPP_MODULE(stan_fit4%model_name%_mod){ &rstan::stan_fit<%model_name%_namespace::%model_name%, boost::random::ecuyer1988>::param_oi_tidx) .method("grad_log_prob", &rstan::stan_fit<%model_name%_namespace::%model_name%, boost::random::ecuyer1988>::grad_log_prob) + .method("hessian_times_vector_log_prob", + &rstan::stan_fit<%model_name%_namespace::%model_name%, boost::random::ecuyer1988>::hessian_times_vector_log_prob) + .method("hessian_log_prob", + &rstan::stan_fit<%model_name%_namespace::%model_name%, boost::random::ecuyer1988>::hessian_log_prob) .method("log_prob", &rstan::stan_fit<%model_name%_namespace::%model_name%, boost::random::ecuyer1988>::log_prob) .method("unconstrain_pars", diff --git a/rstan/rstan/inst/include/rstan/stan_fit.hpp b/rstan/rstan/inst/include/rstan/stan_fit.hpp index da830e176..2dd937de8 100644 --- a/rstan/rstan/inst/include/rstan/stan_fit.hpp +++ b/rstan/rstan/inst/include/rstan/stan_fit.hpp @@ -25,6 +25,11 @@ // void R_CheckUserInterrupt(void); +// TODO: where are these includes actually supposed to go? +#include +#include +#include + // REF: cmdstan: src/cmdstan/command.hpp #include #include @@ -1001,7 +1006,7 @@ class stan_fit { get_all_flatnames(names_oi_, dims_oi_, fnames_oi_, true); // get_all_indices_col2row(dims_, midx_for_col2row); } - + stan_fit(SEXP data, SEXP seed, SEXP cxxf) : data_(data), model_(data_, Rcpp::as(seed), &rstan::io::rcout), @@ -1177,6 +1182,113 @@ class stan_fit { END_RCPP } + + /** + * Expose the hessian_log_prob of the model to stan_fit so R user + * can call this function. + * + * @param upar The real parameters on the unconstrained + * space. + * @param jacobian_adjust_transform TRUE/FALSE, whether + * we add the term due to the transform from constrained + * space to unconstrained space implicitly done in Stan. + */ + SEXP hessian_log_prob(SEXP upar, SEXP jacobian_adjust_transform) { + BEGIN_RCPP + std::vector std_par_r = Rcpp::as >(upar); + Eigen::Matrix par_r = + Eigen::Map>( + std_par_r.data(), std_par_r.size()); + + if (par_r.size() != model_.num_params_r()) { + std::stringstream msg; + msg << "Number of unconstrained parameters does not match " + "that of the model (" + << par_r.size() << " vs " + << model_.num_params_r() + << ")."; + throw std::domain_error(msg.str()); + } + //std::vector par_i(model_.num_params_i(), 0); + Eigen::Matrix grad_f; + Eigen::Matrix hess_f; + + double lp; + if (Rcpp::as(jacobian_adjust_transform)) + stan::model::log_prob_hessian( + model_, par_r, lp, grad_f, hess_f, &rstan::io::rcout); + else + stan::model::log_prob_hessian( + model_, par_r, lp, grad_f, hess_f, &rstan::io::rcout); + + Rcpp::NumericVector grad = Rcpp::wrap(grad_f); + Rcpp::NumericMatrix hess = Rcpp::wrap(hess_f); + + hess.attr("log_prob") = lp; + hess.attr("gradient") = grad; + SEXP __sexp_result; + PROTECT(__sexp_result = Rcpp::wrap(hess)); + UNPROTECT(1); + return __sexp_result; + END_RCPP + } + + /** + * Expose the hessian_log_prob of the model to stan_fit so R user + * can call this function. + * + * @param upar The real parameters on the unconstrained + * space. + * @param jacobian_adjust_transform TRUE/FALSE, whether + * we add the term due to the transform from constrained + * space to unconstrained space implicitly done in Stan. + */ + SEXP hessian_times_vector_log_prob( + SEXP upar, SEXP v, SEXP jacobian_adjust_transform) { + + BEGIN_RCPP + + // Is there a way to convert directly from SEXP to an Eigen matrix? + std::vector std_par_r = Rcpp::as >(upar); + Eigen::Matrix par_r = + Eigen::Map>( + std_par_r.data(), std_par_r.size()); + + std::vector v_r = Rcpp::as >(v); + Eigen::Matrix v_vec = + Eigen::Map>( + v_r.data(), v_r.size()); + + if (par_r.size() != model_.num_params_r()) { + std::stringstream msg; + msg << "Number of unconstrained parameters does not match " + "that of the model (" + << par_r.size() << " vs " + << model_.num_params_r() + << ")."; + throw std::domain_error(msg.str()); + } + //std::vector par_i(model_.num_params_i(), 0); + Eigen::Matrix hess_f_dot_v; + + double lp; + if (Rcpp::as(jacobian_adjust_transform)) + stan::model::log_prob_hessian_times_vector( + model_, par_r, v_vec, lp, hess_f_dot_v, &rstan::io::rcout); + else + stan::model::log_prob_hessian_times_vector( + model_, par_r, v_vec, lp, hess_f_dot_v, &rstan::io::rcout); + Rcpp::NumericVector hess_f_dot_v_r = Rcpp::wrap(hess_f_dot_v); + + // TODO: do this differently + hess_f_dot_v_r.attr("log_prob") = lp; + SEXP __sexp_result; + PROTECT(__sexp_result = Rcpp::wrap(hess_f_dot_v_r)); + UNPROTECT(1); + return __sexp_result; + END_RCPP + } + /** * Return the number of unconstrained parameters */