We propose additive regression models for response densities, which may be observed only via a sample from the conditional density given scalar covariates - or even a single observation, nicely linking our functional regression approach to distributional regression models for scalar data. We formulate our models in Bayes Hilbert spaces to incorporate nonnegativity and integration to one constraints for densities. Considering densities and Bayes Hilbert spaces with respect to an arbitrary finite measure allows us to unify models for continuous, discrete (compositional data) and mixed densities.We propose a (penalized) maximum likelihood estimation approach and derive asymptotic results for the estimators, also giving us (asymptotic) inference. We show a computationally useful approximate equivalence of the maximization problem to (penalized) estimation for certain multinomial - or equivalently, Poisson - additive regression models, greatly simplifying our implementation in an R package. We apply our new approach to a motivating question and data set from gender economics, where we model the mixed density of the woman’s share in a couple’s total labor income - which is continuous on (0, 1) with positive probability mass on 0 and 1 for single-earner couples - in Germany, given age of the youngest child, East vs. West German residence, and year.