Dear Pystan Helpers,
From the pystan doc:
The model is defined as:
model_code =
data {
int<lower=0> N;
vector[N] x;
int<lower=0,upper=1> y[N];
}
parameters {
real alpha;
real beta;
}
model {
y ~ bernoulli_logit(alpha + beta * x);
}
We would like to use pystan vb.
sm = StanModel(model_code)
fit_vb = sm.vb(data=data)

We assume that vb uses “meanfield” as default. Just to confirm, does it mean that the variational distribution is a product of normal distributions? How does it work in Pystan? Does it look at the parameters and then introduce a normal distribution for each parameter? Thus, above we would have two normal distributions: one for alpha and one for beta. Is this correct?

How do we extract the parameters for the variational distribution, after vb finishes running?
It looks like we can only get the posteriors samples. We would like to know the parameters as well. 
How can we extract the ELBO value in each iteration?
Thank you.