Random truncated normal variates.

def rtruncated_normal(mu, tau, a=-np.inf, b=np.inf, size=None):
    """
    Random truncated normal variates.
    """

    sigma = 1. / np.sqrt(tau)
    na = utils.normcdf((a - mu) / sigma)
    nb = utils.normcdf((b - mu) / sigma)

    # Use the inverse CDF generation method.
    U = np.random.mtrand.uniform(size=size)
    q = U * nb + (1 - U) * na
    R = utils.invcdf(q)

    # Unnormalize
    return R * sigma + mu