Efficient code to sum two dependent random variables

Thank you! This is all super helpful. Good to learn that one can pass a vector as constraints.

Unfortunately, I’m still a bit stuck. I would like to find the factors that contribute to the onset time. In particular, in my full model in which t_O is integrated out, the variance of the onset distribution is constant (1 parameter, homoscedastic), but the mean of the onset distribution is defined as a function O_{mean}= \alpha + \beta \cdot x where the parameters are the intercept \alpha and \beta coefficients and x is the “design matrix” of observations associated with the predictor variables. So, I am unclear how to constrain the t_O values when they are defined in terms of a function of predictor variables. This is a bit like a simplex in which the set of values is constrained so that the weighted sum of them is less than t. But in this case, I fear the manifold may be a bit more complex than a simplex, but I really haven’t figured this out. I alluded to this issue in my first post, although my thinking has progressed a lot since that time, but the core issue persists. Is there a way to put “joint” constraints on parameters so that a function of the parameters is constrained to be less than a particular value? Or, is there a better way to implement this?