Call MvNormal.random() method inside model

Hmm, your use case does sound reasonable. For censored distribution usually we integrate the truncated tail area, but since it is difficult for multivariate distribution I think your solution is good to try.