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

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.

I don’t believe the variational distribution is exposed anywhere, if that’s what you’re looking for.

It is possible to compute the log_density and gradients of the likelihood itself (https://pystan.readthedocs.io/en/latest/api.html#pystan.StanFit4model.log_prob).