This function generates m sets of imputed values for the outcome of interest under a fitted joint causal model. The
algorithm draws parameters from the posterior distribution of the model which are then used to obtain simulated
responses (from the posterior predictive distribution of the missing values) via a rejection algorithm. The
bound for acceptance/rejection is obtained via a trust region optimisation.
The imputed values are used to create m complete imputed datasets and perform complete data
analysis and inference about the parameters of interest using any summary statistics.