The only feature suggestion I can suggest is the addition of a convolution operator to combine distributions (reified as RVar's in this implementation, though of course the difference between a random variable over a distribution and the distribution is rather thin)
I don't think I understand. My familiarity with probability theory is fairly light. Are you referring to the fact that the PDF of the sum of random variables is the convolution of their PDFs? If so, the sum of random variables can already be computed as "liftA2 (+) :: Num a => RVar a -> RVar a -> RVar a" since RVar is an applicative functor (or using liftM2 since it's also a monad). Or perhaps you mean an operator that would take, say, 2 values of the 'Uniform' data type and return an instance of the 'Triangular' type corresponding to the convolution of the distributions?