I found out this by implementing simple normal inverse function that each is fairly slow to use. It seems that one should avoid each as much as possible to take leverage of the vectorisation in Q.
However, for an example below, taking a Poisson distribution
(exp neg lambda) * prd (x#lambda) % 1 + til x
til itself is not vectorisable… i.e. If i put a vector after til, it will not work.
actually I tried to implement it today. This augmented solution is a beast. The only concern I have is that essentially we are trading off between speed and memory.