How to get the joint distribution in Pyro?

Thanks for the clarifications! But why does model correspond to p(x,z;θ)? Does Pyro do something like producing a histogram from the data x and the latent variable z generated from q(z|x;ϕ) and then using the histogram to approximate p(x,z;θ)? How is z related (linked) to x? After all, to get the ELBO (the loss for optimization), we need approximations for both p(x,z;θ) and q(z|x;ϕ).