Output of log_prob and optimizing not consistent for inverse Wishart model.

Summary:

When obtaining a provisional value for the posterior mode via postmode = optimizing(...), I believe postmode$value should be the same as log_prob(fit, unconstrain_pars(fit, postmode$par), adjust_transform=FALSE). With an inverse Wishart in the model, this is not true. In fact, the difference isn't even constant in the parameter.

Description:

The help file for log_prob states

log_prob 'signature(object = "stanfit")'
Compute the log posterior ('lp__') for the model represented by a
'stanfit' object. Note that, by default, 'log_prob' returns the
log posterior in the unconstrained space; set 'adjust_transform =
FALSE' to make the values match Stan's output.

and the help file for optimizing states

Value:

If the optimization is done successfully, a list with named
components:

par: The point estimate found. Its form (vector or list) is
     determined by the 'as_vector' argument.

value: The value of the log-posterior (up to an additive constant,
     the '"lp__"' in Stan) corresponding to 'par'.

So optimizing(...)$value and log_prob(...) should match if log_prob() is called on optimizing(...)$par and care is taken to transform the parameter to the unconstrained space and not include the Jacobian from the transformation. This does not happen.

In a much more complicated example than below, I created the logposterior function myself, and it was the same as optimizing()$value except for an additive constant, so I think the problem lies in log_prob() or unconstrain_pars(). log_prob seems to work fine on models that do not include an inverse Wishart component, so I suspect that the constraints for a covariance matrix aren't being handled by unconstrain_pars() correctly.

Reproducible Steps:

Stan model code:

data {
  int<lower=0> nobs;
  int<lower=0> ndim;
  matrix[nobs, ndim] y;
  real mu_mn;
  real<lower=0> mu_sd;
  real<lower=ndim> sig_df;
  cov_matrix[ndim] sig_scale;
}
parameters {
  vector[ndim] mu;
  cov_matrix[ndim] sigma;
}
model {
  for(i in 1:nobs){
    y[i] ~ multi_normal(mu, sigma);
  }
  mu ~ normal(mu_mn, mu_sd);
  sigma ~ inv_wishart(sig_df, sig_scale);
}

R code

library(rstan)

set.seed(234234)

nobs <- 1000
ndim <- 10
sigma <- crossprod(matrix(rnorm(2*ndim^2),ncol = ndim))
cholsig <- chol(sigma)
mu <- rnorm(ndim)
y <- matrix(0, nobs, ndim)
for(i in 1:nobs){
  y[i,] <- mu + crossprod(cholsig, rnorm(ndim))
}

standat <- list(y = y, mu_mn = 0, mu_sd = 10, sig_df= ndim + 1, sig_scale= diag(1, ndim))

model <- stan_model("invwish.stan")

fit <- stan("invwish.stan", data = standat, chains = 1, iter = 1)

postmode1 <- optimizing(model, data = standat, as_vector = FALSE)

## iter = 10 to get a different parameter
postmode2 <- optimizing(model, data = standat, as_vector = FALSE, iter = 10)


lp1 <- log_prob(fit, unconstrain_pars(fit, postmode1$par), adjust_transform = FALSE)
lp2 <- log_prob(fit, unconstrain_pars(fit, postmode2$par), adjust_transform = FALSE)


lprobs <- cbind(c(lp1, postmode1$value, lp1 - postmode1$value), c(lp2, postmode2$value, lp2 - postmode2$value))
rownames(lprobs) <- c("log_prob", "optimizing", "diff")
colnames(lprobs) <- c("par1", "par2")

lprobs

Current Output:

 lprobs
                 par1       par2
log_prob   -136468.17  -241054.5
optimizing  -18389.39 -2462650.2
diff       -118078.78  2221595.7

Expected Output:

 lprobs
           par1         par2
log_prob     X            Y
optimizing   X            Y
diff         0            0

or at least

 lprobs
           par1         par2
log_prob     X            Y
optimizing   X + C        Y + C
diff         C            C

RStan Version:

2.14.1

R Version:

R version 3.3.2 (2016-10-31)

Operating System:

Windows 10 Education
Version 1607
OS Build 14393.693

Start by running the supplied inverse-Wishart Stan model and R reproduction, then compare optimizing()$value with log_prob() after unconstrain_pars() using adjust_transform=FALSE. Trace the handling of covariance-matrix constraints and inverse-Wishart evaluation in those entry points; done means the values agree, or differ only by a parameter-independent constant, across the reproduced cases.