Dear pystan helpers,

We have the model:

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)

Is there an easy way to compute log p(z,y) using pystan?

We can compute this manually, but we would like to know if pystan supports this?

How do we define p(z)? We assume that if we use pystan vb, we know the default is meanfield, so we assume p(z) is a product of Gaussians.

Is there an easy way to compute log p(z,y) using pystan?

Thank you.