After having performed a regression, how can I estimate the probability of a specific outcome? Let's assume that I have a dataset which gives me the two variables: height ($x_1$), daily calorie intake ($x_2$) and weight ($y$) of a person.
In this dataset, we assume I have a large enough number ($n$) of examples in my dataset to make it statistically accurate.
Now, I can do a regression with multiple variables (in this case: 2), to obtain a prediction of a person's weight based on their height and calorie intake. It's rather easy: If I want to find out the predicted weight of a person, I need simply insert the value $x_1$ and $x_2$ into my fitted equation to obtain $y$.
However, I am wondering how I can do something a bit different:
Given a height $x_1$ and a daily calorie intake $x_2$, what is the probability of a person having a weight greater than or less than a certain value, or having a weight in a certain range.
For example: Assume I have $2000$ calories intake per day and a body height of $180 \text{ cm}$, what is the probability that the person weighs between $80\text{ kg}$ and $90\text{ kg}$.
Is such an estimation possible? And if so, how can I do it (analytically or numerically)? Is there any way for me to bound the error that this estimate will have?
 A: If you're taking a frequentist approach, it sounds you're trying to compute the distribution that's used in a prediction interval.
You can do that with a regression, but you do rely on the distributional assumption more than for other forms of inference.
The usual way is to construct a version of a pivotal quantity, but there are some differences from the usual definition of pivotal quantity because instead of trying to get a pivot for a parameter, you need something that will work for a future observation.
In regression, if $y_f$ is that future observation and $\hat{y}_f=x_f\hat{\beta}$ its predicted value, then $Q=\frac{y_f-\hat{y}_f}{s_f}$ has a t-distribution with $n-p-1$ d.f. where $p$ is the number of predictors and $s_{f}$ is the standard error of the prediction.
The formula for $\text{Var}(y_f-\hat{y}_f)$ is often given in derivations of the prediction interval; $s_f$ is the square root of its sample estimate. If $x_f$ is the row-vector of predictors (including the constant) for the future observation, $y_f$, then
$\text{Var}(y_f-\hat{y}_f) = \sigma^2(1+x_f (X'X)^{-1} x_f')$
To obtain $s_f$, replace $\sigma^2$ by $s^2$ in the above formula and take the square root.
You can make probability statements about $Q$, and from those, derive statements about $y_f$ in terms of $\hat{y}_f$ and $s_f$ that should have the appropriate coverage if the assumptions hold.
Packages that will give one-sided prediction intervals can be used to help with the sort of computation you need, but it's easy enough to do directly.
A: You can do it easily if you know the distribution of the error terms (or make an assumption). In case you describe, normality could be reasonable. Then $Y \sim N(\beta_1 x_1 + \beta_2 x_2, \sigma^2)$. Lets say that $\beta_1=0.4$,  $\beta_2=0.004$ and $\sigma^2 = 100$. Then for the height and intake you describe, the weight is distributed as $Y \sim N(80, 10^2)$, and hence probability that his weight is between 80-90 kilos is roughly $\Phi(1) - \Phi(0)\approx 0.34$). But be advised that there are multiple assumptions which should be checked so that this inference makes sense, such as constant variance and the normality or errors.
