# A simple multivariate normal distribution failed. Anything wrong?

Hello,
I am testing my toy code for a multivariate normal distribution.
After making a simple toy data, I did inference.
Simply
x\sim \text{MVN}(\mu,\Sigma)
\mu \sim N(3,5)

Also, I put some general priors over \Sigma.
I think this is a very simple problem, but the MCMC fails.
I am wondering if I did any mistake or it is nature of multivariate normal distribution.
I put error messages and seesion info below. But, it looks like all the 4 chains diverge neither they are similar.

Thanks

library(rstan)
library(mvtnorm)

set.seed(1234)
cov=matrix(c(1,0.5,0.3),nrow=3,ncol=1)%*%matrix(c(1,0.5,0.3),nrow=1,ncol=3)
# create data from
dat=mvtnorm::rmvnorm(n=100,mean=c(0,6,7),sigma=cov)

stan_code='data {
int<lower=1> N; // number of observations
int<lower=1> J; // dimension of observations
vector[J] x[N]; // observations
}
parameters {
corr_matrix[J] Omega;
vector<lower=0>[J] sigma;
vector[J] mu; // a vector of Zeros (fixed means of observations)
}
transformed parameters {
cov_matrix[J] Sigma;
}
model {
mu~normal(3,5);
sigma ~ lognormal(0, 3); // prior on the standard deviations
Omega ~ lkj_corr(2); // LKJ prior on the correlation matrix
x ~ multi_normal(mu,Sigma); // sampling distribution of the observations
//target += multi_normal_lpdf(x|mu,Sigma);
}'

dat2=list(x=dat,N=dim(dat),J=dim(dat))
fit2=stan(model_code=stan_code,data=dat2,cores=4,iter = 5000,chains=4)


Warning message (BTW, I didn’t try change tree-depth etc… as it is fairly simple problem).

Warning messages:
1: There were 157 divergent transitions after warmup. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
2: There were 9797 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
3: There were 3 chains where the estimated Bayesian Fraction of Missing Information was low. See
http://mc-stan.org/misc/warnings.html#bfmi-low
4: Examine the pairs() plot to diagnose sampling problems

5: The largest R-hat is NA, indicating chains have not mixed.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#r-hat
6: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#bulk-ess
7: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#tail-ess

Session info

R version 4.0.4 (2021-02-15)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.2 LTS

Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0

locale:
 LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8 LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8
 LC_MESSAGES=C.UTF-8 LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C LC_TELEPHONE=C
 LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C

attached base packages:
 stats graphics grDevices utils datasets methods base

other attached packages:
 bridgesampling_1.0-0 bayesplot_1.8.0 mvtnorm_1.1-1 rstan_2.21.2 ggplot2_3.3.3

loaded via a namespace (and not attached):
 Brobdingnag_1.2-6 tinytex_0.30 tidyselect_1.1.0 xfun_0.22 purrr_0.3.4 lattice_0.20-41
 V8_3.4.0 colorspace_2.0-0 vctrs_0.3.6 generics_0.1.0 stats4_4.0.4 loo_2.4.1
 utf8_1.2.1 rlang_0.4.10 pkgbuild_1.2.0 pillar_1.5.1 glue_1.4.2 withr_2.4.1
 DBI_1.1.1 plyr_1.8.6 matrixStats_0.58.0 lifecycle_1.0.0 stringr_1.4.0 munsell_0.5.0
 gtable_0.3.0 codetools_0.2-18 coda_0.19-4 labeling_0.4.2 fst_0.9.4 inline_0.3.17
 callr_3.5.1 ps_1.6.0 parallel_4.0.4 curl_4.3 fansi_0.4.2 Rcpp_1.0.6
 scales_1.1.1 RcppParallel_5.0.3 jsonlite_1.7.2 farver_2.1.0 gridExtra_2.3 digest_0.6.27
 stringi_1.5.3 processx_3.5.0 dplyr_1.0.5 grid_4.0.4 cli_2.3.1 tools_4.0.4
 magrittr_2.0.1 tibble_3.1.0 crayon_1.4.1 pkgconfig_2.0.3 ellipsis_0.3.1 Matrix_1.3-2
 prettyunits_1.1.1 ggridges_0.5.3 rstudioapi_0.13 assertthat_0.2.1 R6_2.5.0 compiler_4.0.4

Your covariance matrix is of rank one:

> cor(dat)
[,1] [,2] [,3]
[1,]    1    1    1
[2,]    1    1    1
[3,]    1    1    1


In retrospect that is clear from the code, but it was well camouflaged. ;)
(Unrelated to that, an efficient implementation would use Cholesky and an isotropic normal.)

1 Like

And the SUG provides some input also :)

1 Like

Thanks a lot. My stupid mistake. I had a wrong equation for cov matirx. Now I think I got correct inference…

set.seed(1234)
mat_temp=matrix(rnorm(9),3,3)
mat_temp[upper.tri(mat_temp)]=0
cov=mat_temp%*%t(mat_temp)
#cov=matrix(c(1,0.5,0.3),nrow=3,ncol=1)%*%matrix(c(1,0.5,0.3),nrow=1,ncol=3)
dat=rmvnorm(n=100,mean=c(0,6,7),sigma=cov)

1 Like