Allowing for matrix parameter model declarations

[Sliderland]

Posterior objects currently store parameter dimensions using named integers as follows:

dimensions = list("alpha" = 1,
          			 "beta" = 3
                              )

This works for scalar or vector-valued parameters, but fails for matrix parameters. For example, consider the following matrix declaration in stan:

matrix[2, 3] beta

One might expect to be able to store this type of parameter in the posterior object as follows

dimensions = list("alpha" = 1,
          			 "beta" = c(2, 3)
                                )

Which, in fact, the current implementation allows you to do. However, when computing posterior reference draws, the current implementation of posterior_reference_names does not allow for this. It requires that every entry be an integer. In order to abide by this restriction, one could declare each element of the matrix as its own parameter:

dimensions = list("alpha" = 1,
						"beta[1,1]" = 1,
						"beta[1,2]" = 1,
						"beta[1,3]" = 1,
						 ...
						"beta[2,3]" = 1	
					)

But this could be extremely tedious for larger models, especially one with many different matrix parameters. It would also obscure a lot of the useful information this list provides.

Finally, one could imagine that we should just be able to declare beta = 6, but posterior reference draws checks will be expecting beta[1], beta[2], beta[3], beta[4], beta[5], beta[6] to be in the draws list. Stan will instead return beta[1,1], beta[1,2], ..., causing an error after sampling.

The attached file contains the workflow of adding a simple model which contains a matrix parameter. In order to run this model for themselves, one must simply Line 7 to point to their local instance of posteriordb. I would recommend creating a new local branch of posteriordb before running.

The main areas of interest for this file are the comments in Lines 137-152, which are the relevant sections for defining the posterior's model dimensions, and Lines 205-206, where we use posteriordb to compute the reference draws. The dimension list in Lines 154-165 would be correct under the current implementation of posteriordb.

library(rstan)
library(posteriordb)

options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)

Sys.setenv(PDB_PATH = "/Users/Work/Documents/posteriordb/")
pdbl <- pdb_local()

added_by <- "GitHub Issue"
added_date <- Sys.Date()

data_name <- "syn_multi_normal"
model_name <- "multi_normal_unknown_mean_cov"
reference_name <- paste(data_name, model_name, sep = "-")

s <- "
data {
  int<lower=1> N;
  int<lower=1> K;
  array[N] vector[K] y;
}

parameters {
  vector[K] mu;
  cov_matrix[K] Sigma;
}

model {
  mu ~ normal(0, 1);
  Sigma ~ inv_wishart(
    K + 2,
    diag_matrix(rep_vector(1, K))
  );
  y ~ multi_normal(mu, Sigma);
}
"

set.seed(123)

data <- list(
  N = 100L,
  K = 3L,
  y = cbind(
    rnorm(100, mean = 0, sd = 1),
    rnorm(100, mean = 2, sd = 4),
    rnorm(100, mean = 4, sd = 0.5)
  )
)

model <- rstan::stan_model(
  model_code = s,
  model_name = model_name
)

draws <- rstan::sampling(
  model,
  data = data,
  seed = 123
)

# ------------------------------------------------------------------
# Add the data
# ------------------------------------------------------------------

data_info <- as.pdb_data_info(
  list(
    name = data_name,
    title = "Synthetic Three-Dimensional Normal Data",
    description = paste(
      "Representative synthetic data consisting of 100 observations",
      "of three continuous variables."
    ),
    keywords = c("synthetic", "multivariate"),
    urls = NULL,
    references = NULL,
    added_by = added_by,
    added_date = added_date
  )
)

pdb_data <- as.pdb_data(
  data,
  info = data_info
)

write_pdb(
  pdb_data,
  pdbl,
  overwrite = TRUE
)

# ------------------------------------------------------------------
# Add the model
# ------------------------------------------------------------------

model_info <- as.pdb_model_info(
  list(
    name = model_name,
    title = "Multivariate Normal with Unknown Mean and Covariance",
    description = paste(
      "A multivariate normal model with an unknown mean vector",
      "and covariance matrix."
    ),
    keywords = c("multivariate", "covariance_matrix"),
    framework = "stan",
    urls = NULL,
    references = NULL,
    added_by = added_by,
    added_date = added_date
  )
)

pdb_model <- as.model_code(
  model,
  info = model_info
)

write_pdb(
  pdb_model,
  pdbl,
  overwrite = TRUE
)

# ------------------------------------------------------------------
# Add the posterior
# ------------------------------------------------------------------

pdb_posterior <- as.pdb_posterior(
  list(
    pdb_model_code = pdb_model,
    pdb_data = pdb_data,
    keywords = c("multivariate", "covariance_matrix"),
    urls = NULL,
    references = NULL,
    
    # What one might expect to implement:
    #
    # dimensions = list(
    #   mu = 3,
    #   Sigma = c(3, 3)
    # ),
    #
    # OR
    #
    # dimensions = list(
    #   mu = 3,
    #   Sigma = 9
    # ),
    
    # Under the current implementation, how one could implement
    # this type of model
    
    dimensions = list(
      "mu" = 3,
      "Sigma[1,1]" = 1,
      "Sigma[1,2]" = 1,
      "Sigma[1,3]" = 1,
      "Sigma[2,1]" = 1,
      "Sigma[2,2]" = 1,
      "Sigma[2,3]" = 1, 
      "Sigma[3,1]" = 1, 
      "Sigma[3,2]" = 1, 
      "Sigma[3,3]" = 1
    ),
    
    reference_posterior_name = reference_name,
    added_by = added_by,
    added_date = added_date
  )
)

write_pdb(
  pdb_posterior,
  pdbl,
  overwrite = TRUE
)

# ------------------------------------------------------------------
# Add the existing Stan draws as reference draws
# ------------------------------------------------------------------

reference_info <- as.pdb_reference_posterior_info(
  list(
    name = reference_name,
    inference = list(
      method = "stan_sampling",
      method_arguments = list(
        chains = 10,
        iter = 20000,
        warmup = 10000,
        thin = 10,
        seed = 123
      )
    ),
    diagnostics = NULL,
    checks_made = NULL,
    comments = "Posterior meant to show .",
    added_by = added_by,
    added_date = added_date,
    versions = NULL
  )
)

reference_draws <- compute_reference_posterior_draws(reference_info, pdbl)
reference_draws <- check_reference_posterior_draws(reference_draws)

write_pdb(
  reference_draws,
  pdbl,
  overwrite = TRUE
)