mean() over only some distributions, letting others remain as distributions

[rossry]

I'm trying to handle a situation where I'm modeling both population heterogeneity and distribution over uncertainty, something like:

individual_risk = lognormal(0,1) # this is a global variable, but uncertain
susceptibility = lognormal(0,1) # this varies across the population

populationwide_average_risk = mean_but_only_over_susceptibility(individual_risk ^ susceptibility)
# I want this^ to be distributed over uncertainty in individual_risk, but averaged over the population's susceptibility

In some cases, the function will be sufficiently separable that you can take a simple mean() over the part you want to mean over, and leave the other distribution on the outside...but in this case (and I assume, a substantial fraction of cases) that doesn't work, even with a nifty math hack.

I don't care what the syntax is like, nor do I have a great suggestion. If I had to come up with something, I'd probably say something like mean_over([susceptibility], individual_risk ^ susceptibility), where the first argument is a list of the distributions to collapse by meaning.