# Predicting a realistic distribution using regression

I'm trying to predict a distribution of a continuous variable, that looks like the real distributions I see in my training data. As an example, say I'm trying to predict people's wages, and I know various characteristics about each person e.g. age, gender, job type.

I have a large training set, and I have tried to build a model using multiple regression.

Using this model I try to predict wages for new people (i.e. people not in my training data) - from my understanding I am predicting the average wage for a person with those characteristics. So if I plot the distribution of my predictions it does not look my distribution of real wages, as these are averages. How do I go about making my predictions that look more 'real'? It seems to me I have to add in the 'error term', but how do I do this without making assumptions about the errors?

On looking at the errors, they are not normally distributed, so I tried some non-parametric prediction alogorithms e.g. random forest, and decision trees. But the predictions from these models had a similar distribution to the ones from my regression model.

So how do I go about making a realistic distribution from my predictions, without making assumptions about the distribution of the errors?

Sorry if this is a simple question, or I'm misunderstanding - I'm new to machine learning and predictive models. This seems like it should be a common issue, but I can't find any guides that talk about it.

• One possibility is to use tree regression, and then just draw histograms for each of the the data sets defined by the terminal nodes. Fit some curves afterwards if you want to, Mar 5, 2021 at 15:45
• Thanks for the suggestion. I’m not sure I understand - do you know of any examples I could look at?
– rw2
Mar 5, 2021 at 17:56
• It's not hard at all. Each tree defines a subset of the data. Take all the $Y$ values in a given subset and draw a histogram of them. There is one of your density estimates. You will get one histogram for each terminal node. There is R code and examples in Chapter 18 of "Understanding Regression Analysis: A Conditional Distribution Approach." (amazon.com/…) Mar 5, 2021 at 18:20

Yes, predictions from most other prediction methods (like random forests or CARTs) will also only give "averages". More specifically, all these tools output single numbers, or point predictions. They typically aim at giving you the conditional expectation of the outcome, given covariate values.

So what you are looking for is predictive densities. In the time series forecasting community, these are also called density forecasts.

Unfortunately, there are much less theory, work and tools for density forecasting than for (expectation) point forecasting. On the one hand, that is because point forecasting is easier to do. On the other hand, it's because point forecasting is easier to explain.

There are a few ways you can go about this. However, no approach will be completely free of assumptions. Standard Ordinary Least Squares assumes homoskedastic normally distributed errors, so all it needs to estimate is the variance of the error term. (The mean of the error is zero, because we are predicting the expected value.) You have already determined that this does not work for you. An alternative, e.g., for time series forecasting, is autoregressive conditional heteroskedasticity, or ARCH (or GARCH), where the distribution of the errors changes over time.

Another approach would be something like a Poisson or Negative Binomial regression, where the output would be a predicted Poisson or NegBin distribution, and the full distribution will depend on your covariates.

Yet another way is using Neural Networks for density forecasting. There is some work on that, mainly by binning data and then predicting probabilities that an outcome will fall into a given bin.

Finally, you could also fit your model and collect in-sample residuals, then resample from those residuals and add the resamples to the predicted expected values. This assumes that the distribution does not depend on the covariates, but at least it makes no assumptions on the distributional family. You can weaken the homoskedasticity assumption by stratifying the resampling, but then you are of course still making assumptions on the strata.

Which approach makes most sense will depend on what you are actually doing. For instance, retail sales forecasting can be done quite meaningfully using NegBin regression, because observations are overdispersed counts, so this makes intuitive sense. If you have continuous data, then some GLM might be useful, e.g., using a gamma family. You may want to consider asking a more specific question.

• Thanks for the clear and detailed answer, it's been really helpful to read. I'm thinking of adding re-sampled residuals, as you suggest, and maybe grouping those residuals according to co-variates (the dataset is very large, and most covariates are categorical, so hope this will work). I have a question - I'm not sure I understand how Poisson or NegBin regression helps here - is this just GLM with poisson or NegBin family? Does this not also just provide point predictions? Do you know of any links to an explanation?
– rw2
Mar 4, 2021 at 16:19
• Your proposed resampling and grouping by covariate values is exactly the stratified resampling I was referring to, so that makes a lot of sense. Poisson and NB regression fit parameters for a Poisson or NB distribution, and then you can simply draw from this fitted distribution, which will be different depending on the covariate values. I did density forecasting for retail sales in this way (Kolassa, 2016, which may be useful). ... Mar 5, 2021 at 7:28
• ... It's really analogous to fitting with OLS, then resampling from the fitted error distribution. Essentially, in OLS you would resample from $N(X\hat{\beta}, \hat{\sigma}^2)$, in (say) Poisson regression, you would resample from $\text{Pois}(\exp(X\hat{\beta}))$, and in NB regression, from $\text{NB}(\exp(X\hat{\beta}),\hat{o})$ (the overdispersion parameter $o$ is usually not linked to the regressors $X$, though it could easily be, but rather assumed fixed). Mar 5, 2021 at 7:31