# Impute a distribution

I am doing an applied project. I am wondering if there is any way to impute a distribution in Stan.

Separating a larger model into several smaller models can be either a computational approximation (e.g., replace a multilevel model with a two stage regression), or even a theoretical necessity (e.g., the calibration model and regression model can be assumed as independent in a measurement-error model, or the ''do" operator in causal inference.).

``````data{
real y;
}
parameter{
real<lower=0>  sigma;
}
model{
y~normal (beta, sigma);
}
``````

Instead of making inference on the exact posterior distribution of beta given y, suppose in the ideal situation I have already fitted another calibration model and know beta is exactly 1.
Then I can impute beta into the data block

``````data {
...
real beta;
}
``````

More realistically, I may only know the distribution of beta from the separately-fitted calibration model, say I know p(beta | calibration model )= N(0,1). And I assume beta is independent of y, such that p(beta | y)=N(0,1). Conceptually, I am imagining something like

``````data {
...
real beta;
beta~normal(0,1);
}
``````

In other words, it imputes the value of beta by a pre-fixed distribution. If beta is discrete, I know how to do that, as I can marginalize out beta by finite mixture. When beta is continuous, I have to calculate the likelihood as an integration. Setting a prior beta~N(0,1) in the model block does not achieve the same effect.

Arguably, this is not the most Bayesian approach, for I should always update beta after I see y. Nevertheless at least as an approximation inference, is there anyway to impute the distribution (in the posterior distribution)?

1 Like

Canâ€™t you still marginalise it out by Monte Carlo (with a lot of `log-sum-exp`s)? Maybe with the new fancy MPI stuff if the number of imputed draws is big?

Integrating the likelihood Monte-Carloly is one solution, but probably not scalable with the dimension of distribution to be imputed especially when it is done manually. The bias-variance trade-off again: a slow Monte-Carlo integration vs a biased point estimation.
And yes, MPI should give a more elegant approach.

More thoughts: I can always mix HMC with Gibbs sampling (via a loop outside Stan), but that seems to miss the main point. When I donâ€™t impute distributions and make exact HMC sampling, the computation cost is d^0.25, which is in fact less than Gibbs. Shouldnâ€™t I do even better if I know a priori the posterior distributions of some parameters?

Iâ€™m not really sure what youâ€™re looking for here. You need to integrate out the parameter and if all youâ€™ve got is a bag of samples, youâ€™ll need to do Monte Carlo or some variant. If youâ€™ve got a continuous function (like `beta ~ normal(0,1)`) there either is or will soon be an `integrate_1d` function that you could use.

Looks like you would like to have â€ścutâ€ť functionality (see. e.g.
Better together? Statistical learning in models made of modules). WinBUGS has this (maybe JAGS, too), but itâ€™s not possible in Stan directly.

You can reduce the required number of Monte Carlo draws by various quadrature rules up to few dimensions, CCD up to 10-15 dimensions, and quasi Monte Carlo for more dimensions.

1 Like