Hi folks,

I have a 2D bimodal posterior distribution and I’m trying to sample from this posterior.

Your help would be gratefully appreciated if you could guide me whether there is any way i could sampling from both modes using Stan.

Here is the model information:

Given two correlated Gaussian random variables, u1 and u2 with zero mean, unit variance and correlation coefficient rho equal to 0.9, consider the following transformation:

x = u1* a

y = (u2/a)+b*(u1^2+a^2)

where a = 1.15 and b=0.5.

I used the standard normal bivariate distribution “z” as my prior.

and i believe I’ve done all of the transformations correctly.

I used 8 different chains for my simulation.

The problem is that in all chains, it is just sampling from one mode at the time and cannot sample from both modes together.

I would gratefully appreciate it if you could guide me in this regard.

I attach the simulation results for 8 chains as well. (The plot axis are “x” and “y” variables and the curve line in plot is the “g” function)

Here is my Stan code:

data{

vector[2] mu;

matrix[2, 2] Sigma;

real a;

real b;

matrix[2, 2] L;

real J;}

parameters {

vector[2] z;

}transformed parameters {

// z is bivariarte standard normal, u is correlated bivariate normal and L is lower triangular Cholesky matrix.

vector[2] u = L*z;

real x = a * u[1];

real y = (u[2]/a) + (b * (u[1]^2+(a^2)));

real g = (10^2) - (0.5 * (x+y)^2) - (((x-y)^2)/0.5);

}model {

//prior (Bivariate standard normal)

z ~ multi_normal(mu, Sigma);

//likelihood

target += normal_lpdf(g | 0,0.15)// J is Jacobian transformation but it has constant value.

target += log(J);

}generated quantities {

int I = 0;

{

real gg = (10^2) - (0.5 * (x+y)^2) - (((x-y)^2)/0.5);

if (gg <= 0) I = 1;

}

}

Thank you so much for your time and consideration.

Hamed