Is it possible to access the Hessian from the output when using cmdstan’s optimize method in R or to specify a draws argument so I can get estimates of standard errors on the constrained parameter space?

I was using the optimize function from rstan, but when I moved over to using my universities’ HPC I ran into issues with rstan. I’ve found a workaround with CmdstanR but I can’t seem to figure out where in the output object the Hessian might be stored, or how to return it, or how to get draws from the constrained parameter space.

Here’s how I would get what I want from rstan::optimizing():

```
library(cmdstanr)
library(rstan)
# Use example model from cmdstan
file <- file.path(cmdstan_path(), "examples", "bernoulli", "bernoulli.stan")
mod <- stan_model(file)
stan.data <- list(N = 10, y = c(0,1,0,0,0,0,0,0,0,1))
# Get optimization results
stan_output <- optimizing(mod,
data = stan.data,
algorithm = "LBFGS",
hessian = TRUE,
draws = 50000)
# Extract parameter names
par_names <- names(stan_output$par)
# Get parameter estimates
par_values <- stan_output$par
# By grabbing many draws of the parameters, we can take the square root
# of the variances in the cov matrix for these draws to estimate SE
par_se <- sqrt(diag(cov(stan_output$theta_tilde)))
# Get our t-values
t_val <- par_values / par_se
output <- data.frame(parameters = par_names,
est = par_values,
se = par_se,
t_val = t_val,
row.names = NULL)
```